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Ph! Abstract 

2 . It was believed until very recently that a near-equatorial satellite would always keep 

^ \ up with the planet's equator (with oscillations in inclination, but without a secular 

■ drift). As explained in Efroimsky and Goldreich (2004), this misconception originated 

from a wrong interpretation of a (mathematically correct) result obtained in terms 
^ , of non-osculating orbital elements. A similar analysis carried out in the language of 

^-^ ' osculating elements will endow the planetary equations with some extra terms caused 

^ by the planet's obliquity change. Some of these terms will be nontrivial, in that 

they will not be amendments to the disturbing function. Due to the extra terms, 
the variations of a planet's obliquity may cause a secular drift of its satellite orbit 
inclination. In this article we set out the analytical formalism for our study of this drift. 
We demonstrate that, in the case of uniform precession, the drift will be extremely slow, 
because the first-order terms responsible for the drift will be short-period and, thus, 
will have vanishing orbital averages (as anticipated 40 years ago by Peter Goldreich), 
while the secular terms will be of the second order only. However, it turns out that 
variations of the planetary precession make the first-order terms secular. For example, 
the planetary nutations will resonate with the satellite's orbital frequency and, thereby, 
may instigate a secular drift. A detailed study of this process will be offered in a 
subsequent publication, while here we work out the required mathematical formalism 
and point out the key aspects of the dynamics. 



* In this article, we use the word "precession" in its most general sense which embraces the entire spectrum 
of changes of the spin-axis orientation ~ from the long-term variations down to the Chandler-type wobble 
down to nutations and to the polar wonder. 
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1 Physical motivation and the statement of purpose 



Ward (1973, 1974) noted that the obhquity of Mars may have suffered large-angle motions 
at long time scales. Later, Laskar and Robutel (1993) and Touma and Wisdom (1994) 
demonstrated that these motions may have been chaotic. This would cause severe climate 
variations and have major consequences for development of life. 

It is a customary assumption that a near-equatorial satellite of an oblate planet would 
always keep up with the planet's equator (with only small oscillations of the orbit inclina- 
tion) provided the obliquity changes are sufficiently slow (Goldreich 1965, Kinoshita 1993). 
As demonstrated in Efroimsky and Goldreich (2004), this belief stems from a calculation 
performed in the language of non-osculating orbital elements. A similar analysis carried 
out in terms of osculating elements will contain hitherto overlooked extra terms entailed by 
the planet's obliquity variations. These terms (emerging already in the first order over the 
precession- caused perturbation) will cause a secular angular drift of the satellite orbit away 
from the planetary equator. 

The existence of Phobos and Deimos, and the ability of Mars to keep them close to 
its equatorial plane during obliquity variations sets constraints on the obliquity variation 
amphtude and rate. Our eventual goal is to establish such constraints. If the sateUites' 
secular inclination drifts are slow enough that the satellites stay close to Mars' equator 
during its obliquity changes through billions of years, then the rigid-planet non-dissipative 
models used by Ward (1973, 1974), Laskar & Robutel (1993), and Touma & Wisdom (1994) 
will get a totally independent confirmation. If the obliquity-variation-caused inclination 
drifts are too fast (fast enough that within a biUion or several biUions of years the satellites 
get driven away from Mars' equatorial plane), then the inelastic dissipation and planetary 
structure must play a larger role than previously assumed. 

Having this big motivation in mind, we restrict the current article to building the required 
mathematical background: we study the obliquity-variation-caused terms in the planetary 
equations, calculate their secular components and point out the resonant coupling emerging 
between a satellite's orbiting frequency and certain frequencies in the planet axis' precession. 
A more thorough investigation of this interaction will be left for our next paper. 

2 Mathematical preliminaries: 

osculating elements vs orbital elements 

Whenever one embarks on integrating a satellite orbit and wants to take into account direc- 
tion variations of the planet's spin, it is most natural to carry out this work in a co-precessing 
coordinate system. This always yields orbital elements defined in the said frame, which are 
ready for immediate physical interpretation by a planet-based observer. A well camoufiaged 
pitfall of this approach is that these orbital elements may come out non-osculating, i.e., 
that the instantaneous ellipses (or hyperbolae) parametrised by these elements will not be 
tangent to the physical trajectory as seen in the said frame of reference. 
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2.1 The osculation condition and alternatives to it 



An instantaneous orbit is parametrised with six independent Keplerian parameters. These 
include the three Eulerian angles i, u, Q which define the orientation of the instantaneous 
orbital plane, the eccentricity and semimajor axis of the instantaneous ellipse or hyperbola, 
and the mean anomaly at an epoch, , which determines the initial position of the body. 
Often are employed sets of other variables. The variables are always six in number and are 
some functions of the Keplerian elements. These are, for example, the Delaunay set, the 
Jacobi set, and two sets offered by Poincare. More exotic is the Hill set which includes the 
true anomaly as one of the elements (Hill 1913). 

Systems of planetary equations in terms of the above sets of parameters may be de- 
rived not only through the variation-of-parameters (VOP) method but also by means of the 
Hamilton- Jacobi one. However, the latter approach, though fine and elegant, lacks the power 
instilled into the direct VOP technique and cannot account for the gauge invariance of the 
N-body problem (Efroimsky 2002a,b), important feature intimately connected with some 
general concepts in ODE (Newman & Efroimsky 2003). The Hamilton- Jacobi technique 
implicitly fixes the gauge, and thus leaves the internal symmetry heavily veiled (Efroimsky 
& Goldreich 2003). Below we shall present several relevant facts and formulae, while a more 
comprehensive treatment of the matter may be looked up in (Efroimsky & Goldreich 2003, 
2004). 

As well known, a solution 

r = /(Ci,...,C6,t) (1) 

to the reduced two-body problem 

•■ Gm r 

^ + ^ - = ^ (2) 

is a Keplerian ellipse or hyperbola parameterised with some set of six independent orbital 
elements Cj which are constants in the absence of disturbances. 

In the N-particle setting (or, more generally, in the presence of whatever disturbances) 
each body will be subject to some perturbing force AF(r, r, t) acting at it from bodies 
other than the primary: 

•J, Gm r ^ -1 

? + ^ - = AF , 3 

Solving the above equation of motion by the VOP method implies the use of ([T]) as an ansatz, 

r = / (Ci(t), G,{t), t) , (4) 

the function / being the same as in ([T]), and the "constants" Gi now being endowed with 
a time dependence of their own. Insertion of (jlj) into (|3]) is insufficient for determining the 
six functions Gi{t) . because the vector equation ([3]) comprises only three scalar equalities. 
To furnish a solution three more equations are needed. Since the age of Lagrange it has been 
advised in the literature to employ for this purpose the conditions of osculation. 



imposition whereof ensures that the physical velocity 

dt dt dG. dt 



— * 

dr df ^ df dGi , ^ 

^ + E 7^ , (6) 
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in the disturbed case is expressed by the same function g(Ci, ... , C^, t) as it used to in the 
two-body configuration: 

^ - f « 

The conditions being essentially arbitrary, their choice affects only the mathematical devel- 
opments but not the physical orbit. The orbit's invariance under the alternative choices of 
the supplementary constraints refiects the gauge freedom, i.e., the internal symmetry present 
in this problem. In the language of pure mathematics this is a fiber-bundle structure. In 
physicists' terms it is an example of the gauge freedom. Essentially, this is merely a case of 
ambiguous parameterisation. 

The Lagrange constraint and the law of motion, Q, with ansatz (jl]) inserted therein, 
yields the following equation of the elements' evolution: 

? 1^" § - 1: • 

[Cn Cj] being the matrix of Lagrange brackets introduced as 

— » — » 

df dg df dg 

So defined, the brackets depend neither on the time evolution of Cj nor on the choice of 

— * — * 

supplementary conditions, but solely on the functional form of / {Ci^,,,^ , t) and g = df /dt. 
In case we decide to relax the Lagrange constraint and to substitute for it 

$ being some arbitrary function of time and parameters Ci (but, for the reasons of sheer 
convenience, not of time derivatives of Cj), then instead of ([8]) we shall get, for n = 1, ... , 6 , 



dCn dCj / dt dCn dCn dt dCr, 



$ . (11) 



These gauge-invariant perturbation equations of celestial mechanics, generated by direct 
application of the VOP method, were derived in (Efroimsky 2002b). They give us an op- 
portunity to use, in an arbitrary gauge the Lagrange brackets ([9]) defined in terms of 

— * 

the unperturbed functions / and g . Expressions for these brackets have long been known, 
and they can be employed to write down the planetary equations in their gauge-invariant 
form. If our goal were to arrive to the customary Lagrange system of planetary equations, we 
would now stick to the trivial gauge (|5]) and restrict ourselves to conservative forces solely: 
AF = dR{r)/dr. We however bear in mind the objective of exploring the interplay of two 
freedoms: freedom of reference frame choice and that of gauge fixing. For this reason not 
only shall we leave the gauge arbitrary but shall also reserve for our disturbance a form 
general enough to include both physical and inertial forces, the latter showing themselves 
en route from an inertial coordinate system to an accelerated one. An example at hand 
is a satellite orbiting a precessing oblate planet: the oblateness will give birth to an extra 
physical force, while the precession will result in inertial forces if we decide to describe the 
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orbit in non-inertial axes associated with the precessing planet. What is important about 
the inertial inputs is that they are not only position but also velocity dependent. This gives 
us a motivation to consider velocity-dependent disturbances. Another motivation for this 
comes from the relativity where the correction to Newton' gravity law contains velocities 
(Brumberg 1992) 

Expression (1111) is the most general form of the equations for the orbital elements, in 
terms of the disturbing force. To get the generic form of the equations in terms of the 
Lagrangian disturbance, let us recall that the (reduced) two-body problem is described by 

the unit-mass Lagrangian £o(?, r, t) = r /2 — U{y , t), canonical momentum p = r , 
and Hamiltonian "Hoi?, P; ^) = P + f/(f , t) , while the disturbed setting will be 
described by the perturbed functions: 

2 

£(r, F, t) = £o + A£ = ^ - f/(r, t) + A£(r, F, t) , (12) 



P = ? + , (13) 

or 



and 



-I 2 

n{r, p, t) = p? - £ = ^ + t/ + AH 
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(14) 



AH(r, p,t) = - AC 



1 (dACV 

2 [~W 



ATi being a variation of the functional form: AT-L = 7i{r, p, t) — Hoi^, P, t) ■ The 
Euler-Lagrange equations written for the perturbed Lagrangian will read: 

dU 

? = - ^ + AF , (15) 
or 

the new term 

dr dt \ dr J ^ ^ 

being the disturbing force. We see that, in the case of velocity-dependent perturbations, this 
force is equal neither to the gradient of the Lagrangian's perturbation nor to the gradient 
of negative Hamiltonian's perturbation. This should be taken into account when comparing 
results obtained by different techniques. For example, in Goldreich (1965) the word "dis- 
turbing function" was used for the negative perturbation of the Hamiltonian. The gradient 
of so defined disturbing function was not equal to the disturbing force. 

Substitution of f|T6|) into ffTTl) then yields the generic form of the equations in terms of 
the Lagrangian disturbance^ The result. 



Direct plugging of (IT6l) into pT|) results in 



^ ^ " dt dCn dr dCn dt[ J dCn 

J 
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Of d (dAC 



+ $ 



d 



dCn 



AC 



1 
2 



'dACV 
, dr J 



dCj 



df d 



[17) 



dCn dCn dt 




dAC' 
dr , 



not only reveals the convenience of the generalised Lagrange gauge 

- dAC 
$ = 



dv 



(which reduces to $ = in the case of velocity- independent perturbations), but also explic- 
itly demonstrates how the Hamiltonian variation comes into play: it is easy to notice that 
the sum in square brackets is equal to — A'H^™"*-'. Substitution of the general- type force 
AF into f|TT]) or of the Lagrangian perturbation AL into f lT7|) gives us the gauge- invariant 
planetary equations, provided we know the expressions for elements of the inverse to the 
Lagrange-bracket matrix [Cj Cj] . When the elements Cj are chosen as the Kepler or 
Delaunay variables, (fT7|) entails the gauge-invariant versions of the Lagrange and Delaunay 
planetary equations (See Appendix 1 to Efroimsky & Goldreich (2003).) 



2.2 Goldreich (1965) 

The earliest attempts to describe satellite motion about the precessing and nutating Earth 
were undertaken by Brouwer (1959), Proskurin & Batrakov (1960), and Kozai (1960). 

In 1965 Goldreich accomplished a ground-breaking work that marked the beginning of 
studies of the Martian satellite dynamics. He started out with two major assertions. One was 
that the Martian satellites had either been formed in the equatorial plane or been brought 
therein very long ago. The second was that Mars has experienced, though its history, a 
uniform precession. While the former proved to be almost certainly correct j§ the validity 
of the latter remains model-dependent. While in the simplest approximation the planetary 



If the disturbance were to depend upon positions solely, the first term on the right-hand side of the above 
equation would be equal simply to 9AL/9C„ . In general, though, we should rely on the chain rule 

9A£ 9A£ df d /^C dr dAC df dAC d{i + 4) 



dCn dr dCn dr dCn dr dCn dr 5C„ 
insertion whereof into the preceding formula entails: 



dt dCn d¥ dCn \ dCn dt dCn / V d'v 



To arrive herefrom to p7)) . one will have to add and subtract [1 /2)d{d{AC) / dr) / dCn on the right-hand 
side. Thereby one makes the gauge function $ appear everywhere in the company of d{AC) / dr . 

^ Phobos and Deimos give every appearance of being captured asteroids of the carbonaceous chondritic 
type (Veverka 1977, Pang et al. 1978, Pollack et al. 1979, Tolson et al. 1978). Given the considerations 
summarised in detail by Murison (1988), it seems least unlikely that the Martian satellites were formerly 
asteroids that were captured by the mechanism of gas drag. If so, this must have occurred early in the 
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precession is uniform, a more involved analysis, carried out by Ward (1973, 1974), Laskar 
& Robutel (1993), and Touma & Wisdom (1994), offers evidence of strongly nonuniform, 
perhaps even chaotic, variations of the Martian obliquity at long time scales. It should 
be said, though, that the analysis presented by these authors was model-dependent. In 
particular, it was performed in the approximation of the planet being a nondissipative rigid 
rotatorjl Since the orbits of both satellites are located within less that 2° from the equator, 
the two said assertions come to contradiction, unless there exists a mechanism constraining 
satellite orbits within the vicinity of the primary's equator. (Otherwise, as Goldreich noted 
in his paper, "the present low inclinations of these satellites' orbits would amount to an 
unbelievable coincidence.") In quest of such a mechanism, Goldreich (1965) investigated 
evolution of the Kepler elements of a satellite in a reference system co-precessing with the 
planet. He followed the traditional variation-of-parameters (VOP) scheme, i.e., assumed a 
two-body setting as an undisturbed problem and then treated the inertial forces, emerging 
in the co-precessing frame, as perturbation (along with another perturbation caused by the 
oblateness of the planet). Below we present a brief summary of Goldreich's results, with 
only minor comments. 

Goldreich began by applying formulae ([12]- [IS]) to motion in a coordinate system attached 
to the planet's centre of mass and precessing (but not spinning) with the planet. In this 

history of the solar system while the gas disk was substantial enough to effect a capture. At that stage of 
planetary formation, the spin of the forming planet would be closely aligned with the orbital plane of the 
planet's motion about the Sun i.e., the obliquity would be small and the gas disk would be nearly coplanar 
with the planetary orbit. Energetically, a capture would be most likely to be equatorial. This is most 
easily seen in the context of the restricted three-body problem. The surfaces of zero velocity constrain any 
reasonable capture to occur from directions near the inner and outer collinear Lagrange points (Szebehely 
1967, Murison 1988), which lie in the equatorial plane. Also, a somewhat inclined capture would quickly be 
equatorialised by the gas disk. If the capture inclination is too high, the orbital energy is then too high to 
allow a long enough temporary capture, and the object would hence not encounter enough drag over a long 
enough time to effect a permanent capture (Murison 1988). Thus, Phobos and Deimos were likely (in as 
much as we can even use that term) to have been captured into near-equatorial orbits. 

In all these studies the planet was modelled as a rigid body. Touma & Wisdom (1994) emphasised 
in their research that they assumed Mars to be always in its principal rotation state. A common feature 
of all these models was their neglect of the inelastic relaxation of the spin axis' precession. Whether the 
latter simplification is physically justified is yet to be determined. Inelastic dissipation is, essentially, the 
internal friction and results from alternating stresses and strains emerging in a wobbling rotator (Prendergast 
1958; Burns & Safronov 1973; Efroimsky 2001, 2002c; Molina, Moreno & Martinez-Lopez 2003; Sharma, 
Burns & Hui 2004). It is well-known that the lower the frequency of precession the slower dissipation of the 
elastic energy associated therewith. This circumstance suggests that slow obliquity changes will result in 
virtually no damping. (The same circumstance makes one think that the fastest modes of spin precession 
are well dissipated. It was for this reason that Touma & Wisdom (1994) accepted the model of Mars always 
being in its principal-rotation state.) The real physical picture is far more complicated, mainly due to the 
essentially nonlinear nature of the phenomenon. Since Mars is not a perfectly symmetrical oblate spheroid 
but is a triaxial body, its free rotation in a spin state, which is even slightly different from the principal 
one, is described by the elliptic functions of Jacobi whose decomposition over sines and cosines contains 
an infinite multitude of harmonics. From here, it is possible to demonstrate that in the course of rotation 
a continuous exchange of energy among these harmonics is taking place. Due to this phenomenon, the 
alternating stresses associated with short-period precession are always getting energy from stresses caused 
by long-period changes of the spin axis. This way, dissipation of the short-time motions indirectly leads 
to damping of the long-period ones. This phenomenon, called "energy cascade" is very well known in the 
theory of turbulence, but in fact it is generic and shows itself in nonlinear situations. In the precessing-top 
context it has not been explored so far, and we do not know how effective it is in restraining the obliquity 
variations caused by solar torque. 
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system, the equation of motion includes inertial forces and, therefore, reads: 

/Xxr — /2x(/Xxr) = 



(19) 



dUo d AU 



2/7xr — /Jxr — /Tx(/2xr) 



dr dr 

dots denoting time derivatives in the co-precessing frame and /x standing for the co- 
ordinate system angular velocity relative to an inertial framejll Here the physical (i.e., 
not associated with inertial forces) potential U (r) consists of the (reduced) two-body part 
Uo{r) = —GMf/r^ and a term AU{r) caused by the planet's oblateness. The overall 
disturbing force on the right-hand side of the above equation is generated, according to f|T6|l . 
by 



AC (r, r, = - AU (r) + r ■ (/X x r) + i (/x x r) ■ (/X x r) 



(20) 



Since in this case 



then 



dAC 



dr 



/X X r 



r + 



dAC 
dr 



/X X r 



(21) 



and, therefore, the appropriate Hamiltonian variation will look: 



An 



AC 



1 (dAC' 
2 



dr 



(22) 



AL/ + p ■ (/x X r) ] = AU - (r X p) ■ /x 



AU 



i ■ fl . 



This way, Al-i becomes expressed through quantities defined in the co-precessing frame: 
the satellite's orbital momentum vector J = r x p and the precession rate p, . 

Goldreich employed the above expression in the role of a disturbing function R in the 
planetary equations: 

(U _ cos? d ( - AH) 1 d ( - AU) 

dt na?{l — e^y/'^ sin ? du na?{l — e^)^/^ sin ? dVt ' 



dn ^ 1 d{-AH) ^ ^24) 

dt na?{l — e^y/'^ sin « di ' 



* Be mindful that fl , though being a precession rate relative to an inertial frame, is a vector defined in 
the co-precessing frame. (For details see section 8.6 in Marsden and Ratiu (2003) or section 27 in Arnold 
(1989).) In this frame, 

_ _ din . dhp , „ dhp 



ip and hp being the inclination and the longitude of the node of the planetary equator of date relative to 
that of epoch. 
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where 

— ATi = R = Roblate + Rinertial (25) 

consists, according to fl22|) . of two inputs]^ 

. G Tfl J2 r 2 2/ \ 1 

Rohiate[i^) = -A(7 = — - — — 1 - 3 sin 2 sin [u + u) , (26) 



and 

Rinertial = 3 pL = \J G m O (1 - 6^) • /X . (27) 

Here m = {mprimary + '^secondary) ■ The mean motion is, as ever, n = (G my^'^ a^^^'^ , 
while p stands for the mean radius of the primary, u denotes the true anomaly, and 



e2 



1 + e cosu 



(28) 



is the instantaneous orbital radius. In the right-hand side of ( |27|) it was assumed that the 
angular momentum is connected with the orbital elements through the well-known formula 



J = r X p = yCm a (1 — e^) w (29) 

where 

w = Sti sin i sin — X2 sin i cos i7 + X3 cos i 

is a unit vector normal to the instantaneous ellipse, expressed through unit vectors xi, X2; ^3 
associated with the co-precessing frame Xi, Xg, % (the axes Xi and lying in the planet's 
equatorial plane). The afore- written expression for w evidently yields: 



Rinertial = y G m a (1 — e^) ( /ii sini sini7 — /i2 sini cosfi -|- /is cosi ) , (30) 

while f l26|) may be, in the first approximation, substituted with its secular part, i.e., with its 
average over the orbital period: 

J2 2 3 cos^i - 1 

(^ofeZate) - P _ ^2^3/2 , (31) 

the averaging having been carried out through the medium of formula ( 1109^ from the Ap- 
pendix, with ( l28l) inserted. With the aid of ( l30l) and ( 13T|) . the planetary equations ( HM and 
f l25l) will simplify to: 

di 

— = — /ii cos Q — p2 sinQ , (32) 



^ Our formula ([26]) slightly differs from the one employed by Goldreich (f965), because here we use the 
modern definition of J2 '■ 



r 



m=2 



a being the satellite's latitude in the planet-related coordinate system. The coefScient J used by Goldreich 
(f965) differs from our J2 by a constant factor: J = (3/2) J2 p^/r^ . 
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The latter results in the well-known node-precession formula, 

^] = fi. - ^nJ2 (^)' ^^^^ (t-Q . (34) 

Its insertion into the former entails 

z = - ^ cos [ - X (t - to) + ^^o] + — sin [ - X (t - to) + ^o] + io , (35) 

where 

3 f P^\'^ cos % 

In Goldreich (1965), equation ( l35l) was the main result, its derivation being valid for wobble 
which is slow <^ t/2 ^ ) and close to uniform ( <^ J2 n). 

Despite a warning issued by Goldreich in his paper, this result has often been misinter- 
preted and, therefore, misused in publications devoted to satellites and rings of wobbling 
planets, as well as in the literature on orbits about tumbling galaxies. 

In (l35|l i stands for the inclination defined in co-precessing axes associated with the 
planet's equator, and therefore ( l35l) clearly demonstrates that, in the course of obliquity 
changes, this inclination oscillates about zero, with no secular shift accumulated. Does this 
necessarily mean that the satellite orbit, too, oscillates about the equatorial plane, without 
a secular deviation therefrom? Most surprisingly, the answer to this question is negative. 
The reason for this is that so calculated orbital elements, though defined in the co-precessing 
frame, are not osculating therein. In other words, in the frame where the elements are 
introduced, the instantaneous ellipses parameterised by these elements are not tangent to 
the physical orbit as seen in this frame. 

This circumstance was emphasised by Goldreich, who noticed that formula (!29|) normally 
(i.e., when employed in an inertial frame) connects the osculating elements defined in that 
frame with the angular momentum fxp defined in the same, inertial, frame (i.e., with 
fx r). Since in the above calculation the frame is not inertial (and, therefore, the angular 
momentum is different from f x f but is equal to fxp = fx{f'+p,xr)), then the 
orbital elements returned by fl29|) cannot be osculating in this framej§ On these grounds 
Goldreich warned the reader of the peculiar nature of the elements used in his integration. 

To this we would add that it is not at all evident that the inertial-forces-caused alteration 
of the planetary equations should be achieved through amending the disturbing function with 

^ Were these elements osculating in the frame wherein they had been defined, then formula would 
read: r x f — y^Gm a (1 — e^) w , i.e., would connect the elements with the velocity in that frame. In 
reality, though, it reads: r x p — ^ Gm a (1 — e^) w , i.e., connects the elements with the momentum 
p = r + fl X r which happens to coincide with the satellite's velocity relative to the inertial axes. This 
situation was formulated by Goldreich in the following terms: the orbital elements emerging in the above 
derivation are defined in the co-precessing frame, but are osculating in the inertial one. This illustrative 
metaphor should not, though, be overplayed: the fact that the elements emerging in Goldreich's computation 
return the inertial-frame-related velocity does not mean that this inclination may be interpreted as that 
relative to the invariable plane. (The elements were introduced in the co-precessing frame!) 
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the momentum-dependent variation of the negative Hamiltonian, 3 ■ p, . While the common 
fallacy identifies the disturbing function with the negative Hamiltonian perturbation, in real- 
ity this rule-of-thumb works (and yields elements that are osculating) only for disturbances 
dependent solely upon positions, not upon velocities (i.e., for Hamiltonian perturbations 
dependent only upon coordinates, but not upon momenta). Ours is not that case and, 
therefore, more alterations in the planetary equations are needed to account for the frame 
precession, if we wish these equations to render osculating elements. However, if one neglects 
this circumstance and simply amends the disturbing function with J ■ /x , then the planetary 
equations will give some elements different from the osculating ones. It will then become an 
interesting question as to whether such elements will or will not coincide with those rendered 
by ( l28l) when this formula is used in non-inertial frames. 

All these subtle issues get untangled in the framework of the gauge formalism. Appli- 
cation of this formalism to motions in non-inertial frames of references was presented in 
Efroimsky & Goldreich (2004). The main results proven there are the following. 

— * 

1. If one attempts to account for the inertial forces by simply adding the term J ■ p, 
to the disturbing function, with no other alterations made in the planetary equations, then 
these equations indeed do render quantities that may be interpreted as some orbital elements 
(i.e., as parameters of some instantaneous conies). These elements are NOT osculating and, 
therefore, the instantaneous conies parameterised by there elements are not tangent to the 
physical orbit. Hence, these elements cannot, generally, be attributed a direct physical inter- 
pretation0 except in situations where their deviation from the osculating elements remains 
sufficiently small. 

2. By a remarkable coincidence, these non-osculating elements turned out to be identical 
with those emerging in formula f l29p . This coincidence was implicitly taken for granted by 
Goldreich (1965), which reveals his truly incredible scientific intuition. 

3. To build up a system of planetary equations that render osculating elements of the 
orbit as seen in the co-precessing coordinate system, one has not only to add 3 ■ ft to the dis- 
turbing function, but also to amend each of these equations with several extra terms. Some 
of those terms are of order ( \ jl\/{J2 n) )^ ; some others are of order |/2|/(|/x| J2 n) . Most 
importantly, some terms are of first order in the precession-caused perturbation |/2 1/( J2 n) , 
which means right away that the non-osculating elements used by Goldreich (1965) differ 
from the osculating ones already in the first order. While a more comprehensive account 
on this topic, with the resulting equations, will be offered at the end of this section, here 
we shall touch upon only one question which gets immediately raised by the presence of 
such first-order differences. This question is: what are the averages of these differences? 
Stated alternatively, do the secular components of the said non-osculating elements differ 
considerably from those of the osculating ones? Goldreich (1965) stated, without proof, that 
the secular components differ only in high orders over the velocity-dependent part of the 
perturbation. In our paper we shall probe the limits for this assumption. 



''When the orbit evolution is sufBciently slow, the observer can attribute some physical meaning to 
elements of the osculating conic. For example, whenever an observer talks about the inclination or the 
eccentricity of a perturbed orbit, he naturally implies those of the osculating ellipse or hyperbola. 
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2.3 Brumberg and Kinoshita 



A development, part of which was similar to that of Goldreich (1965), was independently 
carried out by Brumberg, Evdokimova & Kochina (1971), who studied orbits of artificial 
lunar satellites in a coordinate system co-precessing with the Moon. In that article, too, the 
non-osculating nature of the resulting orbital variables did not go unnoticed. The authors 
called these variables "contact elements" and stated (though never proved) that these vari- 
ables return not the correct value of the velocity but that of the momentum. Later, one of 
these authors rightly noted in his book (Brumberg 1992) that the contact elements differ 
from the osculating ones already in the first order over the velocity-dependent part of the 
perturbation. In subsection 1.1.3 of that book, he unsuccessfully tried to derive analytical 
transformations interconnecting these sets of variables]^ 

A similar attempt was undertaken in a very interesting article by Ashby & Allison (1993). 
Though the authors succeeded in many other points, their attempt to derive formulae for 
such a gauge transformation was not successfuljj 

The setting, considered by Goldreich (1965) in the context of Martian satellites and by 
Brumberg et al (1971) in the context of circumlunar orbits, later emerged in the article by 
Kinoshita (1993), who addressed the satellites of Uranus. 

Kinoshita's treatment of the problem was based on the following mathematical construc- 
tion. Denote satellite positions and velocities in the inertial and in the co-precessing axes 
with {r ', v'} and with { r , v} , correspondingly!^ Interconnection between them will 
be implemented by an orthogonal matrix A , 

r = Ar' , v = f=*=ir' + if=*' = ii~^r + iv' = - /txr + ip' , (37) 

ft being the precession rate as seen in the co-precessing coordinate system, and the inertial 
velocity v ' being identical to the inertial momentum p ' . Kinoshita suggested interpreting 
this interconnection as a canonical transformation between variables {'f^ ' , p '} and 
{r,p}, implemented by generating function 

F2 = p- Ar' = (A^p) ■ r' . (38) 

* Contrary to the author's statement, formula (1.1.41) m Brumberg (1992) is not rigorous, but is vahd 
only to first order. (To make it rigorous, one should substitute everywhere, except in the denominators, r 
with r — dR/df .) Besides, the author did not demonstrate his derivation of formula (1.1.43), for Mo , 
from (1.1.42). (In Brumberg's book the mean anomaly is denoted with I , not with M . ) 

Most importantly, the qualitative reasoning presented by the author in the paragraph preceding formula 
(1.1.43) is unrigorous and essentially incorrect. The cause of this is that the author compares the planetary 
equations for contact elements, written in a precessing frame, with the equations for osculating ones, written 
in an inertial frame, instead of comparing two such systems (for contact and for osculating elements), both 
of which are written in a precessing frame. This makes a big difference because, as we already explained 
above, transition to a precessing frame does not simply mean addition of an extra term to the disturbing 
function. 

Despite all these mathematical irregularities, the averaged system of planetary equations (1.1.44), "de- 
rived" by Brumberg for the first-order secular perturbations, turns out to be correct in the limit of uniform 
precession. Just as in the preceding subsection we had a reason to praise the unusual intuition of Goldreich, 
so here we have to pay tribute to the excellent intuition of Brumberg, intuition which superseded his flawed 
mathematics. 

^ To carry out the gauge transformation, the authors used a set of intermediate variables {Q^o) ' (o) J" ' 
which were canonical and, at the same time, osculating. As follows from the theorem proven by Efroimsky 
& Goldreich (2003), these variables are nonexistent when the perturbation depends upon velocities. 

We use notations opposite to those in Kinoshita (1993), in order to conform with the notations of 
Goldreich (1965). 
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This choice of generating function rightly yields 



dp 

while the interconnection between momenta will look: 



Af' , (39) 



P' = ^ = A'P. i.e., P=[A') p' = Ap' , (40) 
whence p = \ -\- p, x r . The Hamiltonian in precessing axes will read 

H(r, p) = H^i^'^-'Ur', p') + ^ = (f^', p" ) + ^ . Ar' 

at 

= Hi^nert) (^'^ _ ^ . (fl x V) = H'^''"''''^ {v ' , p ') - {v X p) ■ . (41) 

The Hamiltonian perturbation, caused by the inertial forces, is — (rxp) ■/x=— Jx/T, 
vector J being the orbital angular momentum as seen in the co-precessing frame. Comparing 
this with ([22]), we see that employment of the above canonical transformation is but another 
method of stepping on the same rake. In distinction from Goldreich (1965) and Brumberg 
et al (1971), Kinoshita in his paper did not notice that he was working with non-osculating 
elements. 

The problem with Kinoshita's treatment is that the condition of canonicity in some sit- 
uations comes into contradiction with the osculation condition. In other words, canonicity 
sometimes implicitly contains a constraint that sometimes is different from the Lagrange con- 
straint ([5]). This issue was comprehensively elucidated in the work Efroimsky & Goldreich 
(2003). The authors began with the reduced two-body setting and thoroughly re-examined 
the Hamilton- Jacobi procedure, which leads one from the spherical coordinates and the cor- 
responding canonical momenta to the set of Delaunay variables. While in the undisturbed 
two-body case this procedure yields the Delaunay variables which are trivially osculating 
(and parameterise a fixed Keplerian ellipse or hyperbola), in the perturbed case the situa- 
tion becomes more involved. According to the theorem proven in that paper, the resulting 
Delaunay elements are osculating (and parameterise a conic tangent to the perturbed tra- 
jectory) if the Hamiltonian perturbation depends solely upon positions, not upon momenta 
(or, the same: if the Lagrangian perturbation depends upon positions but not upon veloci- 
ties). Otherwise, the Delaunay elements turn out to be non-osculating (and parametrise the 
physical trajectory with a sequence of non-tangent conies). As one can see from the above 
equation, the Hamiltonian perturbation, caused by the inertial forces, depends upon the 
momentum, and this circumstance indicates the problem. This trap, in which many have 
fallen, is of special importance in General Relativity, because the relativistic corrections to 
the equations of motion are velocity-dependent 



In an interesting article (Chernoivan & Mamaev 1999), the authors addressed the two-body problem on 
a curved background. The curvature entailed a velocity-dependent relativist correction, which was treated 
as a perturbation. After carrying out the Hamilton- Jacobi development, the authors arrived at canonical 
variables analogous to the Delaunay elements. Orbit integration in terms of these variables would be as 
correct as in terms of any others. The problems began when the authors used these elements to come to 
some conclusions regarding the perihelion precession. Those conclusions need to be reconsidered, because 
they were rendered on the basis of Delaunay elements that were non-osculating. Similar comments may be 
made about the work by Richardson & Kelly (1988) who addressed, using the Hamiltonian formalism, the 
relativist two-body problem in the post-Newtonian approximation. 
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3 Planetary equations 



In this section we shall briefly spell out some results obtained in Efroimsky & Goldreich 
(2004) and shall use these results to derive the Lagrange- type planetary equations fl5^- [M|) 
for osculating elements in a coordinate system co-precessing with an oblate primary. 



3.1 Planetary equations for contact elements 

Above, in subsection 2.1, we provided a very short account of the gauge formalism. Ex- 
pression ( !T7|) . presented there, is the most general form of the planetary equations for an 
arbitrary set of six independent orbital elements, written in terms of an arbitrary disturbance 
of the Lagrangian. 

When the elements Ci are chosen to be the Keplerian or Delaunay sets of variables, 
we arrive at the gauge-invariant versions of the Lagrange or Delaunay planetary equations, 
correspondingly. They are written down in the Appendix to Efroimsky & Goldreich (2003). 
Interplay between the gauge freedom and the freedom of frame choice is explained at length 
in Section 3 of the article Efroimsky & Goldreich (2004), which addresses orbits about a 
precessing planet. It is demonstrated in that work that, if one chooses to describe the 
motion in terms of the non-osculating elements that were introduced in a co-precessing frame 
and defined in the generalised Lagrange gaugj^ f lTS]) . then the corresponding Hamiltonian 
perturbation will read: 

A-H(-"*) = - [i?,,;,,,(/) + /x-(/xg)] , (42) 
while the planetary equations ( !T7|) acquire the form 

dC- d ( - AH(™"*) 

or: 

[Cr ^ = ^ [ RoMateif) + M " (/x g) ] , (43) 

— * 

where / and g stand for the undisturbed (two-body) functional expressions of the position 
and velocity via the time and the chosen set of orbital elements: 

r = /(Ci,...,C6,t) 

(44) 

V = g(Ci,...,C6,t) = ^m,...,C6,t) 

— ♦ 

(so that f {Ci(t), ...jCeit), t) and g{Ci(t), ...,0^(1), t) become the ansatz for solving the 
disturbed problem). For RoUate in P2l) - (143!) . one can employ, dependent upon the desired 
degree of rigour, either the exact expression or its orbital average (PTl) . 
In the generalised Lagrange gauge fll8l) the canonical momentum becomes: 

p = r + = g + $ + = g , (45) 

or or 



It is an absolutely crucial point that choice of a gauge and choice of a reference frame are two totally 
independent procedures. In each frame one has an opportunity to choose among an infinite variety of gauges. 
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which means that its functional dependence upon the time and the chosen set of orbital 
elements is the same as it used to be in the unperturbed case where both the velocity and 
the momentum were simply equal to g{Ci, ...,Cq, t) . This also means that A'H*^™"*'' 
coincides with Goldreich's A'H given by (!26|) . 

We see that the generalised Lagrange gauge (ITSl) singles out the same set of non-osculating 
elements which showed up in the studies by Goldreich (1965) and Brumberg, Evdokimova 
& Kochina (1971), - the set of "contact elements." This is why in f l42|) the Hamiltonian 
perturbation was written with the superscript "cont." 

In P3|) the Lagrange-bracket matrix is defined in the unperturbed two-body fashion 
and can, therefore, be trivially inverted. Hence, when the elements are chosen as the 
Keplerian ones, the appropriate equations look hke the customary Lagrange-type equations 
(i.e., like and (125|) above), with the disturbing function given by (l26ll or, equivalently, 

by m- 

da ^ ^ c}(-A?^(-*)) 

dt ~ na dMo ' ^ ^ 



de l-e^ <9(- A7{(^°"*)) (i _ e2)i/2 9 ( - AH^"""*) 



dt na? e dMo na^ e du 



(47) 



dcj - cos Z d(- /\H^-ont)\ ,^ _ ^2^1/2 d ( - /\H^ 



(cont) 



dt na'^{l — e^)^/^ sin z di na^e de 



dt _ cost 9(-AH(-"*)) 1 9(-AH(-"*) 



(48) 



dt na^(l — e^)^/^ sin z dou na'^{l — e^)-*^/^ sin i dCl 



(49) 



dn 1 d(- AH^^""*) 



dt na'^{l — e^)^/^ sin z di 



(50) 



dM, _ 1 - <9(-A7{(-"*)) 2 d [- AH 



(cont) 



dt na^ e de na da 



(51) 
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3.2 Planetary equations for osculating elements 



When one introduces elements in the precessing frame and also demands that they osculate 

— * 

in this frame (i.e., makes them obey the Lagrange gauge $ = 0) then the Hamiltonian 
variation reads0 



while equation (fT7|) becomes: 



Roblatei^^) + M • (/ X g) + (/X X /) ■ (/X X ] 



(52) 



[Cn Ci 



dCi dA'H'-""'^ 



dt 



dC„ 



(53) 



dCr, 



/ X 



dC„ 



d 



To ease the comparison of this equation with f l43|) . it is convenient to split the expression for 
^q^{osc) ^ j^^Q p^^^g. 



Roblateif, t) + fl- if X 



and 



- (/Tx /) ■ (/7x /) , 
and then to group the latter part with the last term on the right-hand side of fl55]) : 



(54) 



(55) 



[Cn Ci 



dt 



dC„ 



(56) 



^ df ^ ^ di 
+ M - I T^TT X g - / X 



One other option is to fully absorb the 0(|/lp) term into A?/ , i.e., to introduce the 
effective "Hamiltonian" 



A^(e//) 



1 - - 

Roblate{y) + P-if y<i) + -iflxf)-{flxf) 



(57) 



and to write down the equations like this: 

Ci] — = + ■ I X g - / X 



dt 



dCr. 



dCn 



dC„ 



A^ - / X 



dCr, 



(58) 



Just as AH^''""') in ([42]), this Hamiltonian variation is still equal to - R{f, t) + fl ■ J = 
— R{f , t) + /7 • (/ X p) . However, the canonical momentum now is different from g and reads as: 
P = g + (m X /) ■ 



16 



For Ci being chosen as the Keplerian elements, inversion of the Lagrange brackets will 
yield the following Lagrange-type system: 



da 2 
dt na 



(59) 



de 1 — 



dt 



na? e 



d ( - AH^^^f^) 



(60) 



(1 - e^y/^ 

na^ e 



d ( - AH^ 



eff) 



duj — cos i 

dt na^{l — e^y/'^ sini 



d ( - A?{(^//)) 



di 



+ 



na^ e 



oe \ oe oe \ oe 



di cos i 

dt na^ (1 — e^y/'^ sin i 



na^ (1 — e^y/"^ sin i 



dj-A-H^^ff^) ^ (df ^ ^ 



/ X 



dl 

on 



dQ _ 1 

dt na^ (1 — e^y/^ sin i 



di 
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dMo 



1 - 
na^ e 



de 



df 

^■1 



de 



de 



(64) 



2 

na 



d ( - An^'ff^ 



da 



oa 



^ da 



terms fl ■ {df /dMo) x g — (dg/dMo) x / ) being omitted in (^9\ - |60]), because these 
terms vanish identically (see the Appendix). In equations fl58l - [64]) . A'H^'^^-^'' is given by 
( 1571) . In a rigorous analysis, with the true- anomaly dependence of all terms in (l59l - [6^ taken 
into account, the Hamiltonian will look, according to (1261-1271) and ( 1301) : 



A^(e//) 



Roblateil^) + ■ (/ X g) + ^ X /) • (m X /) 



Gm J2 
2 ^ 



1 — 3 sin^ i sin^(a; + u) 



G m a (1 - e2) w-pi~]-(p,^f)-(p,^f) 



G m J2 / 1 + e cos z/ 



2 a3 



1 - e2 



1 — 3 sin^ z sin^ (u + u) 



(65) 



— \ Gm a {1 — e^ 



/ii sini sinf2 — /X2 sini cosf2 + /X3 cosi 



1 
2 



X /) ■ (/^ X /) 



Otherwise, when all terms on the right-hand side of ( IS^ - IMI) are substituted with their 
orbital averages (denoted with the (...) symbol), the Hamiltonian will be, due to ( 13T|) : 



(i?oWate) + M ■ (/ X g) + - ( (/X X /) ■ (/7 X /) ) 



(66) 



G m J2 3 cos^ z — 1 
4 a3 (1 _ e2^3/2 



Gma (1 — e^) ( /ii sini sinf2 — /i2 sinz cosfi + p^ cosi 



^ {{fix f). {fix f)) 



18 



The expression for / x g is true-anomaly-independent and, therefore, does not need to 
be bracketed with the averaging symbols (...). The expression for (a* x ■ (^/X x 
through the orbital elements is too cumbersome, and here we do not write it down explicitly. 
(See formula ( I186P in the Appendix.) When we permit ourselves to neglect the 0(|/2p) 
inputs, all three Hamiltonians, /S.'H'y"^^) , /^-^Cconi) ^ AT-L^^^^^ coincide: 



— * 

{Roblate) + M • (/ X 



GmJ2 p2 3 cos^ i - 1 /- ■ ■ ■ n ■ ■ o , 

— TTTT — a/ G m a 1 1 — I I Ui sm z sm i i — U2 sm 2 cos \l + cos 

4 a3 (1 _ e2)3/2 V V / VA- a- a-^ 

Two important issues should be dwelt upon at this point. 

First, we would remind the reader that the function A'H*^'^""*) , given by expression 
yields the correct functional form of the Hamiltonian only in case we express the 
Hamiltonian through the contact elements (and calculate these through (143|) or (146) - [5T]) ). 
In the currently considered case of osculating elements, this A'H'-™"'*'' is NOT the correct 
expression for the Hamiltonian. The correct functional dependence of the Hamiltonian upon 
the osculating elements, ATi^"^^^ , is given by formula f l52p . Though in this dynamical 
problem the Hamiltonian is unique, its functional dependencies through the contact and 
through the osculating elements differ from one another, one being AH*-™"'*'' as in fH2|) . 
another being A7i^"^^^ as in fl52|) . As for the function Al-L'^^^^^ rendered by fl571) . it is not 
really a Hamiltonian, but is simply a convenient mathematical entity. In the approximation, 
where 0(|/xp) terms are neglected, there is no difference between these three functions. 
Despite this, the 0(|/X|) and 0(|/X|) terms do stay in equations (jSHl-Ell) for osculating 
elements. 

Second, we would comment on our use of expressions fl29] - [30]) in our derivation of fl65] 
- |67]). Above, in subsections 2.2 and 3.1, the use of fl27|) and (l30l) was based on formula 
f l2^ which interconnected J = r*xp=/x (^r*+/xx/^ with contact elements 
a , e , i , and Vt . As demonstrated in Efroimsky & Goldreich (2004), in that frame 
r = g — p, X f . Hence, formula fl29|) interconnects the contact elements with 
J = / X • III the present subsection, we use formula (l29|l in its usual capacity, i.e., to 

— * • — * 

interconnect J = rxr = fxg with osculating elements a , e , i , and Q . 
It may seem confusing that, though in both cases this formula can be written down in the 
same way. 



f X g = \l Gm a {I — e^) w , (68) 
its meaning is so different. The clue to understanding this difference lies in the fact that in 

— * 

one case r = g — p, x f (and, therefore, the elements are contact), while in the other 
case r = g (which makes the elements osculating). For more details see Efroimsky & 
Goldreich (2004). 
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4 Comparison of the orbital calculations, performed in 
terms of the contact elements, with those performed 
in terms of the osculating elements. The simplest 
approximation. 

As explained in Section 2, it follows from equations fl23]- that the initially small incli- 
nation remains so in the course of the oblate primary's precession. Whether this famous 
result may be interpreted as keeping the satellites in the near-equatorial zone of a precessing 
planet will depend upon how well the non-osculating (contact) inclination emerging in (1251 - 
|36|1 approximates the physical, osculating, inclination rendered by (159) - |6l]) . 

4.1 The averaged planetary equations 

Comparing equations ([59] - [M]) for osculating elements with equations (jl6] - [51]) for contact 
elements, we immediately see that they differ already in the first order over the precession 
rate /x and, therefore, the values of the contact elements will differ from those of their 
osculating counterparts in the first order, too. A thorough investigation of this difference 
would demand numerical implementation of both systems and would be extremely time- 
consuming. Meanwhile, we can get some preliminary estimates by asking the following, 
simplified, question: how will the secular, i.e., averaged over an orbital period, components 
of the contact and orbital elements differ from one another? To answer this question, we 
perform the following approximations: 

1^ In equations f[59]- [M]) . we substitute both the RoUate term in the Hamiltonian and 
the /x-dependent terms with their averages (so that, for example, the RoUate term will be 
now substituted with (Robiate) expressed via f[3T]) ). 

2. We neglect the terms of order fi ^ . This way, we restrict the length of time scales 
involved. (Over sufficiently long times even small terms may accumulate to a noticeable 
secular correction.) However, we can now benefit from the approximate equality ( 167)) . 

So truncated and averaged system of Lagrange-type equations will read: 




(69) 




(70) 



na^ e 



df 



] 
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duj 



cos I 



dt no?'{l — e^y/"^ sin i 



+ 



(1 - e^f'^ 
na^ e 



+ (/^-^xg-Zx^ ) - (/J /X ^ ) 



de 



de 



de 



de 



(71) 



di cos i 

dt na^ (1 — e^)-*^/^ sin « 



(72) 



(1 - e2)V2 si 



sm z 



di-An^'^ff^) (df - dS\ . I ^ df , 

+ (/^-iT^Xg-Zx^l ) - (/X /X ) 



91] 



91] 



on 



dVL 



dt na^ (1 — e^yi'^ sin ^ 



, (73) 



1 - 

na^ e 



de 



^ , df ^ - 9g 



(74) 



2 



d ( - AU^^ff^ 



da 



the Hamiltonian here being approximated by (!67j) . and the angular brackets signifying or- 
bital averaging. In equations ( 169|) . (!70|) and ( !72|) we took into account that the averaged 
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Hamiltonian fl67p bears no dependence upon and u . (This, though, will not be the 
case for the exact, i/-dependent, Hamiltonian given by (152!) and (126!) ! ) 

Calculation of ( [df /dCj) x g — (dg/dCj) x / ) and — p, f x dg/dCj ) takes 
pages of algebra. A short synopsis of this calculation is offered in the Appendix. Here follows 
the outcome: 



da 



2 



iGm (1 - e2 



(75) 



dl 

de 



X g - / X 



dg 

de 



na^ (3 e + 2 cos u + e'^ cos 
(1 + e cosu) \/\ — 



(76) 



X g - / X - 1 = -2 /.^ 



n a 



1 + e cosz/ 



(77) 



/ X M 



(78) 



n a 



1 + e COS!/ 



sin i { cos i7 cos [ 2 (w + z/) ] — sin Q cos i sin [ 2 (w + z/) ] } + 



n a? \/l — 
1 + e COS!/ 



e sin i { cosf2 cos (z/ + 2 a;) — 2 sin 17 cos i sin (cj + z/) cosw } 



+ 



yW2 



n a? yjl — 
1 + e cosz/ 



sin z { sin cos [ 2 (cu + z/) ] + cos VL cos z sin [ 2 (a; + z/) ] } + 



1 + e cosz/ 



e sin z { sinQ cos (z/ + 2 a;) + 2 cosfi cosz sin (a; + z/) cos a; } 



+ 



1 + e cos z/ 



sin^ i sin [ 2 (cj + u) 



na^ a/1 — 
1 + e cos z/ 



2 e { - 



sin z/ + sin i cos w sin (ci; 



+ ^) } 
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(79) 



1 + e cos 



— I ( — sin O cos i cos 2uj — cos O sin 2uj ) cos^ u — sin^ ^ ) "I" 



2 sin u cos 1/ ( sin fl cos z sin 2u! — cos f2 cos 2u) )] + 



sin f2 cos z cos 2oj — cos Q sin 2uj ) cos z/ + ^ sin cos i sin 2aj — 2 cos Cl cos^ cj^ sin u | 



/ia V -1- — e r r , . • r-> • r, \ / 2 • 2 \ 

+ U2 i ( COS i I COS z COS 2a; — smU sm 2a; j cos — sm i/ ) + 

1 + e cos i/ L ' ^ ^ 



2 sin 1/ cos u ( — cos ^2 cos i sin 2a; — sin fl cos 2a; ) ] + 



cos cos z cos 2a; — sin sin 2a; ) cos z/ + ^ — 2 sin Q, cos^ a; + cos sin 2a; cos « ^ sin z/ | 



+ 1^3 



1 + e cos z/ 



I sin i cos 2a; 



cos 2a; I cos^ u — siv? v 



) — 2 sin 2a; 



smz/ cosz/ 



+ 



e [ sin i cos 2a; cos v — sin i sin 2a; sin z/ ] } , 



dS_ 
9M„ 



X g — / X - — = , 



(80) 



- fi± an {t - to) VT 



(81) 
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. 2 (1 - e') . 
1 + e cos ly 



(82) 



2^2 



(1 - e^) 



1^1 + e cosu) 



2 ' 



(83) 



(84) 



:i - e 



2^2 



[1 + e cos i^y 



{ /ii [ cos Q cos(a; + i/) — sin Q sin(a; + i/) cos i ] sin(a; + i/) sin i 



+ /t2 [ sin fl cos{uj + u) + cos fl sin(a; + i/) cos i ] sin(a; + i/) sin i 



+ fi3 cos^(a; + u) + sin^(a; + u) cos^ i j | 



^ , r -9/ 



(85) 



2\2 



(1 - e^) 



(1 + e cosi/)^ 



{ /ii [ — cos sin(a; + i/) — sin cos(a; + i/) cos i ] sm{u! + z/) 



+ /i2 [ — sinQ sin(c(; + i/) + cosQ cos(a; + u) cos i ] sm{u) + i/) 



+ /t3 sin(a; + i/) cos(a' + u) sin i } 



(86) 
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/ii , and yUs being the Cartesian components of the precession rate, as seen in the co- 
precessing frame, and fi± being the component of the precession rate, aimed in the direction 
of the angular momentum; it is given by flll7p . 

It may seem strange that the right-hand side of (!76|) does not vanish in the hmit of 
e — i- . The absurdity of this will be easily redeemed by the fact that this term shows up 
only in the equation for du/dt and, therefore, leads to no physical paradoxes in the limit 
of a circular orbit. However, for finite values of the eccentricity, this term contributes to the 
periapse precession. 

Another seemingly calamitous thing is the divergence in (l82l) . This divergence, however, 
entails no disastrous physical consequences, because the term ( l82l) shows up only in the 
planetary equation for dMo/dt and simply leads to a steady shift of the initial condition 
Mo . 

4.2 The case of a constant precession rate. 

The situation might simplify very considerably if we could also assume that the precession 
rate p, stays constant. Then in equations fl69l - rM|) . we would take ft out of the angular 
brackets and proceed with averaging the expressions (df/dCj) x g — / x (dg/dCj) ^ 

only (while all the terms with p, will now vanish). It is, of course, well known that this is 
physically wrong, because the planetary precession has a continuous spectrum of frequencies 
(some of which are commensurate with the orbital frequency of the satellite) Nevertheless, 
for the sake of argument let us go on with this assumption. 

The averaging of (1751) and (|5U|) is self-evident. The averaging of (1751- 17^ is lengthy and 
is presented in the Appendix. All in all, we get, for constant p, : 



df ^ ^ dg \ . df ^ 9g \ 3 G m {1 — e 



o2\ 



M • ( X g - / X ^ j ) = , C, = e,n,co, I, Mo . 

Since the orbital averages ( 188|) vanish, then e will, along with a , stay constant for as long 
as our approximation remains valid. Besides, no trace of p will be left in the equations for 
Q and i . This means that, in the assumed approximation and under the extra assumption 
of constant p , the afore-quoted analysis (1231- [36]) . offered by Goldreich (1965), will remain 
valid at time scales which are not too long. At longer scales (of order dozens of millions of 
years and longer) one has to take into consideration the back reaction of the short-period 
terms upon the secular ones (Laskar 1990). Beside the latter issue, the problem with this 
approximation is that it ignores both the long-term evolution of the spin axis and the short- 
term nutations. For these reasons, this approximation will not be extendable to long periods 

The case of the Earth rotation and precession is comprehensively reviewed by Eubanks (1993). The 
Martian short-time-scale rotational dynamics is of an equal complexity, even though Mars lacks oceans and 
the coupling of its rotation with the atmospheric motions is weaker than in the case of the Earth (Defraigne 
et al 2003, Van Hoolst et al 2000, Dehant et al 2000). 
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of Mars' evolution. This puts forward the bigger question: what maximal amplitude of 
obliquity variations could Mars afford, to keep both its satellites so close to its equatorial 
plane? 

Even in the unphysical case of constant ft the averaged equation (!74|) for the osculating 
Mo differs already in the first order over /x from equation fl48p for the contact Mq . In 
the realistic case of time-dependent precession, the averages of terms containing pi do not 
vanish (except for /x ■ {df/dMo) x g — / x (dg/dMo) ) which is identically nil). These 
terms show up in all equations (except in that for a ) and influence the motion. They will 
be the key to our understanding the long-term satellite dynamics, including the secular drift 
of the orbit plane, caused by the precession jl. 

5 An outline of a more accurate analysis. Resonances 
between the planetary nutation and the satellite or- 
bital frequency. 

Precession of any planet contains in itself a continuous spectrum of circular frequencies 
involved!^ 

/oo 
/7(s) e-''' ds , (89) 
-oo 

some modes being more prominent than the others. For our present purposes, it will be 
advantageous to express the precession rate as function of the satellite's true anomaly: 

/oo 
fliW) e-'^" dW , (90) 
-oo 

W being the circular "frequency" related to the true anomaly u . Needless to say, 
p{t) ftiv) , ii{s) , and fOW) are four different functions. However, we take the lib- 
erty of using the same notation /!(...) because the argument will always reveal which 
of the four functions we imply0 As ever, it is understood that the physical content is 
attributed to the real parts of fi{t) and /x(z^) . 

If we now plug the real part of (pi into (|75] - [HZ]) and carry out averaging in accordance 
with formula (11091) of the Appendix, we shall see that the secular parts of these /T-dependent 
terms do not vanish. They reveal the influence of the planetary precession upon the satellite 
orbital motion. Especially interesting are the resonant contributions provided by the integer 
TV's, because these nutation modes are commensurate with the orbital motion of the satellite. 

For plugging ftiv) into the terms flTSl- ISBl). it will be convenient to rewrite fpi]) as 

m-u) = r \p,^'\W) smiWv) + fL^''\W) cosiWv)] dW . (91) 

JO L J 



A more honest analysis should take into account also the direct dependence of the planet's precession 
rate upon the instantaneous position(s) of its satellite(s): /2 = fl(t ; a , e , , a; , i , Mq) ■ However, the 
back-reaction of the satellites upon the primary is known to be an effect of a higher order of smallness (Laskar 
2004), at least in the case of Mars; and therefore we shall omit this circumstance by simply assuming that 
ft = fl{t) = fl{t{u)). 

Evidently, fl{v) is a short notation for fl{t{v) ) . It is also possible to demonstrate that, in the limit 
of vanishing eccentricity, P-{W) « n fl{W n) and W « s/n . 
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Insertion of the latter into fl75] - [5B]l will still result in extremely sophisticated integrals. For 
example, the term (175ll will have the following orbital average: 





/ - I - o 1 \ ^ /Cm (1 - e2) 




(92) 

^"rfz/ '"^^^^^ sin(iyz/) + /iP(iy) cos(W^i/) 
(1 + e cosz/)^ 

(the averaging rule fllOQp being employed). Evaluation of the two integrals emerging in this 
expression can, in principle, be carried out in term of the hypergeometric functions, but the 
outcome will be hard to work with and hard to interpret physically. Even worse integrals 
will show up in the averages of ([76] - [86]) . 

To get an idea as to how much the terms f[75] - [5^ contribute to the secular drift of 
the satellite orbits from the planet's equator, it may be good to first calculate these term's 
averages under the assumption of small eccentricity. Then, for example, ([92]) will simplify 
to 

(93) 

3 iGm / „n2 /-oo r"^^ 



Air 



(l - J^dW J^" du { fi^l\w) sm{Wu) ( 1 - 2 e cos i/ + 3 cos^ z/ - 



+ /iP(Vr) cos{Wiy) (l - 2e cos z/ + 3 cos^ u + 

We immediately see that the above expression consists of two parts: the non-resonant one 
and the resonant one. 

This topic will be addressed in our subsequent paper, where we shall try to determine 
how much time is needed for these subtle effects to accumulate enough to cause a substantial 
secular drift of the orbit plane. 



6 Conclusions 

In this article we have prepared an analytical launching pad for the research of long-term 
evolution of orbits about a precessing oblate primary. This paper is the first in a series and is 
technical, so we deliberately avoided making whatever quantitative estimates, leaving those 
for the next part of our project. 

The pivotal question emerging in the context of this research is whether the orbital planes 
of near-equatorial satellites will drift away from the planetary equator in the cause of the 
planet's obliquity changes. Several facts have been established in this regard. 

First, the planetary equations for osculating elements of the satellite do contain terms 
responsible for such a drift. These terms contain inputs of the first order and of the second 
order in p, , and of the first order in ft , where fl is the precession rate of the primary. 
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Second, the first-order (but not the second-order) terms average out in the case of a 
constant precession rate, which means that in this case their effect will accumulate only over 
extremely long time scales. We would remind that the sort-period terms of the planetary 
equations do exert back-reaction upon the secular ones. While in the artificial-satellite sci- 
ence, which deals with short interval of time, the short-period terms often may be omitted, in 
long-term astronomical computations (dozens of million years and higher), the accumulated 
influence of short-period terms must be taken into account. A simple explanation of how 
this should be done is offered in section 2 of Laskar (1990). Besides, the contribution of the 
secular second-order terms will, too, accumulate over very large time intervals. 

Third, the first-order drift terms do not average out in the case of variable precession. 
Under these circumstances they become secular. This means, for example, that the turbulent 
history of Mars' obliquity - history which includes both long-term changes (Ward 1973, 1974, 
1979; Laskar & Robutel 1993; Touma & Wisdom 1994) and short-term nutations (Dehant 
et al 2000; Van Hoolst et al 2000; Defraigne et al 2004) - might have lead to a secular 
drift of the initially near-equatorial satellites. If that were the case, then the current, still 
near-equatorial, location of Phobos and Deimos may lead to restrictions upon the rate and 
amplitude of the Martian obliquity variations. To render a judgment on this topic, one 
should compute how quickly this drift accumulates. Very likely, this quest will demand 
heavy-duty numerics. 
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Appendix. 

The inertial terms emerging in the planetary equations 



In this Appendix we write down, as functions of the true anomaly, the p,- and /x- 
dependent terms in fl59l - [6^ . We also calculate their orbital averages in the case of a 
constant precession rate ft. 



A.l. The basic formulae 

Formulae (!55l - [60!) contain the two-body unperturbed expressions for the position and 
velocity as functions of time and the six Keplerian elements, ([1]) and ([7]). To find their 
explicit form, one can employ an auxiliary set of dextral perifocal coordinates q , with an 
origin at the gravitating centre, and with the first two axes located in the plane of the orbit: 

q = { r cos I/, rsini/, 0}'^ = a | cos z/ , sin z/ , | ""^ . (94) 

1 + e cos z/ 



• \ n a sin z/ n a ie + cos z/) 



The corresponding velocities will read: 

n a sin v 

{ — sin z/ , (e + cos z/) , } 



T 



^ O r • X ^ 1 T 



(95) 



VI - e 



2 



the radius in being a function of the semiaxis major a , the eccentricity e , and the 
true anomaly v : 

1 - 

(96) 



1 + e cosz/ 



the true anomaly z/ itself being a function of a , e , and of the mean anomaly M = 
Mo + jl^n dt , where n = (Gm)-*^/^ a~^/^ . Then, in the two-body setting, the position 
and velocity, related to some fiducial inertial frame, will appear as: 

r = f{Ci,...,Ce, t) = R{Q, z, u) q(a, e. Mo , t) , 

(97) 

r = g(Ci,...,C6, t) = R(fi, I, u) q(a, e, , t) , 

R(i7, i, w) being the matrix of rotation from the orbital-plane-related axes q to some fixed 
Cartesian axes (xi, X2, xs) in the said fiducial inertial frame wherein the vectors r and 
r are defined. The rotation is parameterised by the three Euler angles: inclination, i ; the 
longitude of the node, Q ; and the argument of the pericentre, u . Thence, as is well known 
(see, for example, Morbidelh (2002), subsection 1.2), 
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COS Q cos a; — sin n sin u cos i — cos Q sin u — sin Q cos u cos ^ sin fl sin ? 

R = sin f2 cos u + cos fl sin a; cos i — sin sin + cos fl cos u cos z — cos Q sin « 

sin u sin i cos w sin i cos z 

insertion whereof, together with ( IMl) . into the first equation of ( 1971) yields 

1 - e2 

fi = a [ cos cos (u} + iy)— sin f2 sin (a; + u) cos z ] , (99) 

1 + e COS!/ 

1 - 

/2 = a [sinf2 cos (w + u) + cosf2 sin (w + u) cosi] , (100) 

1 + e cosu 

1 - 

/a = a sin (w + u) sin z . (101) 

1 + e cosz/ 

Similarly, substitution of f l^ and f pSjl into the second equation of fj^ entails: 

gi = = [— cosQ sm{u + z/) — sin $7 cos(a; + u) cosi + 

yl — 



e ( — cosQ sinu — sinQ cosu cosi) 



Ti a 

g2 = I - [ — sinf2 sin(ci; + i^) + cosfi cos(ci; + v) cosz + 
vl — 



e ( — sinf2 sinw + cosfi cosw coszl 



(102) 



(103) 



^3 = . = sin z [ cos(a; + z/) + e cosw] , (104) 
vl — 

the subscripts 1,2,3 denoting the x^-, xz components in the fiducial inertial frame 
wherein (1971) is written. 



A. 2. The averaging rule 

From equations 



e + cosz/ Vl - e2 sinz/ 

cosi^ = , smi^ = (105) 

1 + e cosz/ 1 + e cosz/ 



it follows that 



dE ^/\ 



e 



2 



9i/ 1 + e cosz/ 
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(106) 



From the first of formulae f llUSI) and from the Kepler equation one can derive: 

dM 1 - 



dE 1 + e cosi^ 
Together, ffT06|) and ffTOTll entail: 

dM dM dE (1 - e^f^ 



(107) 



whence 



dv dE dv (1 + e cosvf 



dM 



:i08) 



211 Jo 2 7T Jo (1 + e cosz/ 



\2 



Calculation of the integral shows that the right-hand side of the above equation is equal 
to unity, which means that the secular parts should be calculated through the following 
averaging rule: 

^1 _ e^f^ r^- du 



{■■■) ^ ' I •••77- ^ (109) 

2 TT JO (1 + e cosz/) 

In what follows below, we try hard to squeeze the calculations as much as possible, but at 
the same time to leave them verifiable for the interested reader. 



A. 3. Calculation of mItt-xs-Zxt^ 

da da 



As evident from fl94l) and from the first equation of fl97l) . 
da [da Uz/I V^aA,,,M„ « du [daj^^^^^^ a du [daj^^^^^^ 



and, therefore. 



df 1 - 

7fxg = -/xg . (Ill) 
da a 



Similarly, from f l95|) and the second equation of fl971) it ensues that 



di ^ /9g\ ^ /di\ fdi/ 
da \da I \du \da 



t, e, Mo 



(112) 



g (9g dt ( du\ g / Gm -*\ dt [ du\ 
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wherefrom 

/x^ = -;^/xg. (113) 
oa 2 a 

In the undisturbed two-body problem, / x g is the angular momentum (per unit of the 



reduced mass) and is equal to \jGm a (1 — e^) w , where the unit vector 

10 = xi sin z sinfi — X2 sin « cosVt + X3 cos? (114) 

is normal to the instantaneous osculating ellipse, the unit vectors xi, X2, X3 making the 
basis of the co-precessing coordinate system x^, xg, X3 (the axes and belonging to 
the planet's equatorial plane). 

Together, ffTTTj) and ffTT3D will give: 



df ^ r dg 3 - ^ 3 /- — ^ 

Trxg-/x — = — / xg = — JGm a 1 - w 

da oa 2a 2a ^ 



;ii5) 



2 V a 
and, thereby. 



3 /Gm (1 - e2) . _ . 

[ Xi sm I sm i Z — X2 sm i cos i z + X3 cos i 



^ • 1^ 5^ " ^ - ^ " 5^ j = 2 V a ^'''^ 

where 

= /ii sin i sini7 — fi2 sin z cosfi + /is cos z . (117) 

Since, for constant /X , 0116^ is //-independent, then in the uniform-precession case it will 
coincide with its orbital average. 



A A. Calculation of At(^xg-/x^ 

oe de 



Just as in the preceding subsection, one can write: 



:il8) 



2 e + cos V + cos ^ ^ _^ ^ dt /du^ 
(1 + e cosu) (1 — e^) du \de ^ 



t, a, Mo 
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whence 



df 



X 



2 e + cos v + e cos u -* 



de (1 + e cosv) (1 — e'^j 

With dnS]) taken into account, the derivative of the two-body velocity will take the form 

t, a, Mo 



(119) 



de \ de 



du I \de 



1 — ^ dt dv \de] 



t, a, Mo 



(120) 



e ^ / G m -\ dt f du^ 



i, a, Mo 



where 



b = 



121) 



no, X 
= { — cos fl sin a; — sin fl cos u cos z , — sin sin a; + cos fl cos a; cos i , sin i cos a; } 



From this we easily get, with aid of flll4p : 



/ X b 



n a 



1 + e cosz/ 



{ sin Q sin i cos ly , — cos Q sin i cos ly , cos i cos u } 



(122) 



1 + e COS!/ 



cos { sin z cos f2 , — sin i cos f2 , cos i } = w 



T _ n a? a/1 — 



1 + e cos z/ 



cos z/ . 



Together, ffTT9|) . f lT20|) . and f lT22|) yield: 



5/ 

(9e 



2 e + cos z/ + e cos z/ 



'1 + e cosz/) (1 — e^) 1 — 



/ X g - / X b 



^ 3 e + cos z/ + 2 cos z/ /— — _^ n a? \/l — e^ 

— w — ^ A/Gma (1 — e"') — w — cosz/ 



(1 + e cosz/) (1 — e^) 



1 + e cos z/ 



^ (3 e + 2 cos v + e^ cos z/) 
w - 



'1 + e cosz/) a/I — 



whence 



(123) 



df ^ dg 

de de 



/i_L 



no? (3 e + 2 cos z/ + cos z/) 
(1 + e cosz/) Vl — 



With aid of integrals 



2tt 



dv 



2 + e' 



(1 + e cosz/) 



TT 



'1 - e 



2^5/2 



(125) 
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and 

r^'^ cos u du — 3 e 



(1 + e cosuf ^ (1 - e2)^/^ 



(126) 



it is easy to demonstrate that, under the assumption of constant jl , the orbital average of 
ffmj) is nil. 



A. 5. Calculation of At(?^xg-/x^ 

oco oco 



A straightforward calculation, based on formulae (|94l) . ( I95l) . (1971) . and (1981) . gives: 



— xg = — - — = — . e smiZ smi smz/ , (127) 

ou ooo oco 1 + e cosz/ 



2 



f^/* ^ \ c^/s c^/i Vl — e ^ , . 

— xg = — g^ - — g^ = — — ■ e COS 12 sm z smz/ , (128) 

ooo / ao; ou 1 + e cosz/ 



5/ .\ df, dh na'W^ . . 

— xg =— 5,2-— 5(1 = — e cos« smz/ . (129) 

ou j ou ou 1 + e cosz/ 
/ 3 

The two-body angular momentum cannot depend on the argument of the pericentre, 

d_ 
du 



/ X g) = , (130) 



whereby 



— — X g — / X — = — 2 e smiZ smz smz/ , (132) 

ou ou j 1 + e cosz/ 



9/ ^ 9g \ n a -\/l — 

7— X g — J X — = 2 e cosiZ sin z sinz/ , (133) 

ou Ou j 1 + e cosz/ 



^xg-/x— =-2 ^ e cosz smz/ . (134) 

ou ou j 1 + e cosz/ 
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Altogether, these formulae may be written down as: 



Of ^ 7^ dg \ _^ n a? Vl — , , 

7^xg-/x— = - 2 w -— e sinz/ , (135) 

ou ouj I 1 + e cosz/ 

w being the unit vector pointing in the direction of the angular momentum. It is given by 
(1TT4|) . Finally, 



— xg-/x— =-2/ix ^ e smi/ . (136) 

^ OUJ ou I 1 + e cost/ 

Since sinz/ is odd and (1 + e cosz/)~^ is even, it is obvious that, for a constant p, , the 
orbital average of (11361) will vanish. 



A. 6. Calculation of /x(T;7rxg-/x-^ 



Once again, formulae dH]), ([95]), (|97D, and ([98]) yield: 

Of .\ 9/2 9/3 



(137) 



2 ^/f^TI 



n a V i — e 



1 + e cosz/ 



[ cos n cos(a; + z/) — sin Q sin(a; + z/) cos ^ ] cos (a; + z/) sin z 



+ e [cosfi cos((X' + z/) — sini7 sin(u; + z/) cos z 1 cosw sinz 

1 + e cos ^ ^ ^ ^ ^ ^ 



9/ .\ 9/3 9/i 

/ 2 

(138) 



[sinf2 cos(ci; + z/) + cosf2 sin(co' + z/) cosz] cos (w + z/) sinz 

1 + e cosz/ 



n a V i — e 



+ e [sini7 cos( lo + z/) + cosfi sin(t(; + z/) cos z 1 coscj sinz 

1 + e cosz/ ^ ^ ^ 



9/ ^ .\ 9/i 9/2 

/ 3 
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1 + e cos V 



{ [ — sin Q cos(a;+ z/) 



— cos Q sin(c(; + z/) cos i] [ — sin sin {ui + cos Q cos i cos(a; + i/)] — 



[ cos Q cos(a; — sin Q sin(a; + v) cos i ] [ — cos Q, sin (a; + i/) — sin Q cos i cos(a; + v)] } 



n Vl - r r • / \ 

H e — smiZ cos cu + z/ — 

1 + e cos i/ ^ ^ ' 



cosQ sin(a; + i^) cos z ] [— sinQ sino^ + cosQ cos i cos a;] 



— [cosQ cos(a; -\- v) — sinJl sin(a; + u) cosi] [— cosQ sin a; — sinQ cos z cosu; ] } 



1 + e cos J 



— I sin(a' + i/) cos (a; + I/) sin^ i + e sino? cos(a; + i/) — sin(a; + i^) cos a; cos^ i j | , 



(140) 



n 



1 + e cosi/ 



[ — sin(a; + v) sin i ] [ — cos VL sin(a; + i/) — sin VL cos(a; + z/) cos i 



+ 



1 + e cosi/ 



e [— sin(a; + i/) sini] [— cosO sin a; — sinO cos a; cos z 



n 0? Vl — 
1 + e cosi/ 



sm z 



sin^ [u + i/) cos O + sin {u + i/) cos {u + i/) sin Q cos z 



+ 



n Vl — 
1 + e cos 



e sin z [ sin (a; + i/) sin a; cos + sin (a; + v) cos a; sin Q cos z ] , 
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dgi „ <9g 



/3 ^ - /: 



1411 



n a 



1 + e cosu 



sin{u + u) sin z [sinf2 sm{u + u) — cosQ cos{u + u) cosi 



+ 



n a 



1 + e cosu 



e sin{u) + u) sin z [sinfi sinw — cosQ cosu cos z 



df . 



(142) 



na^-y/l — 
1 + e cos 



{ [ cos VL cos(a; + u) 



— sinfi sin(a; + z/) cosi] [— cosfi sin (a; + z/) — sinfi cos? cos(a; + z/)] — 



[sinfi cos(a; + z/) + cosfi sin(a; + z/) cos « ] [sin 12 sin (a; + z^) — cos 12 cosi cos(a; + z^)] } 



+ 



n c? a/1 — 
1 + e cos V 



e { [ cos 12 cos(a; + z/) — 



sin 12 sin(a; + v) cosi ] [ — cos 12 sin a; — sin 12 cos z cos a; 



[sin 12 cos(a; + z/) + cos 12 sin(a; + u) cosi] [sin 12 sin a; — cos 12 cosi coscu] } 



1 + e cos u 



simu + u) cos[u + u) sin i + e 



sinu cos{u} + u) + sin{u + u) cos a; cos^ i | 



In ( 11371 - [T38ll and in (11401 - [T4T]) we benefitted from fs being independent of 12 . Now, by 
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subtracting, accordingly, f lMI]) from ffW]) . and ffW]) from f lT^ . and ffTi^j) from ffT^ . we 
arrive, after some intermediate algebra, at the following three expressions: 



sin z { cos Q cos [ 2 (w + z^) ] — sin Q cos 2 sin [ 2 (w + z/) ] } 

1 + e cosz/ 



?1 CL^ ~~ 6^ 

+ e sin z | cosfi cos (i/ + 2 0;) — 2 sinfi cosz siniu + u) cos a; | , 

1 + e COS!/ ^ ^ ^ ^ ^ ^ 



X g - / X § 1 = (144) 



2 



1 + e cos z/ 



sin 2 { sin Q cos [ 2 (cj + z/) ] + cos Q cos z sin [ 2 (cj + z/) ] } 



+ e sin 2 { sinfi cos (z/ + 2 a;) +2 cos 17 cos i sin (cj + z/) cosw } 

1 + e cos z/ 



?^ X g - / X 1^ 1 = (145) 



n a \/l — . n ■ . r / \ 

sm ^ sm 2 w + z/ 

1 + e cosz/ ^ ^ ' 



n a? Vl - r • / x • / \ 2 • 

+ e <^ 2 smu; cos (z/ + tj) — 2 sm (a; + z/) cosw cos z 

1 + e cosz/ 



na^ a/1 — . 2 ■ • rr^ / \n na^A/1 — ^ r . . o • • / \ 1 

sm I sm 2 (w + z/) J + I e \ — sm z/ + sm i cos w sm (w + z/j ^ 

1 + e cos z/ 1 + e cos z/ J 



Hence, the overall expression for /x ■ y[df/dVL) x g — / x (Sg/Sfi) j will be given by 
( !78|l . To average this expression, in the simple case of a constant p. , it will be sufficient 
to separately average the above three expressions (11431 - I145P and then to multiply the 
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averages by /ii , /i2 and fi^ , correspondingly, and to add everything up. To carry out this 
procedure, it will be convenient first to regroup terms in (11431 - [T^5|) and to through out the 
inputs proportional to the odd functions sin u and sin 2i> , because these inputs will vanish 
after averaging via formula fllOQp . Thus we get: 



sin i { cos cos 2ui cos 2i/ — cos sin 2uj sin 2i/ 

1 + e cosz/ 



sinQ cosz sin 20; cos2i/ + sinf2 cosz cos2co' sin2z/ } 



+ e sin i { cos cos 2w cos v — cos sin 2bj sin z/ 



1 + e cosz/ 



sin Vt cos « sin 2a; cos i/ + 2 sin cos i cos^ w sin v | ) 



n vl — e 
1 + e cos z/ 



2 

sin z ( cos cos 2(jj — sin cos « sin 20; ) ( cos 2z/ + e cos v ) ) 



'1 — e^)^^^ 9 I o I cos 2z/ + e cos z/ 

V 1 — e"* sin « ( cosii cos 2a; — sinii cosz sin 2a; ) / rfz/ = 

Wo n + fi cosz/r 



n a 



2 TT -/o (1 + e cosz/) 



^xg-/x^l )= (147) 



n a? a/1 — 
1 + e cosz/ 



sin « { sin Vt cos 2cc; cos 2z/ — sin f2 sin 2uj sin 2z/ 



+ cosf2 cos « sin 20; cos2z/ + cosfi cosz cos 20; sin2z/ } 
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+ e sin i { sin Q cos 2u cos — sin Q sin 2uj sin z/ 

1 + e cosu 



+ cosfl cosi sin 2a; cosz/ + 2 cosfi cos ^ cos^cj sinz/| ) 



n 0? V^l — 6^ 

sin z (sinfi cos 20; + cosf2 cosi sin 20;) (cos2z/ + e cosi/) ) 

1 + e cos z/ 



I'l — e^)^^^ 9 / ^ / \ /■^'^ cos 2z/ + e cos 1/ , 

V 1 — e"' sm 2 I smiZ cos2o; + cosi 2 cos ^ sin2o; ) / ^ dv = 

Jo (1 + e cosi/r 



n a 



2 IT Jo (1 + e cos u 



(|^xg-/x||| )= (148) 



3 



1 + e cos z/ 



na a/1 — o 

sm i ( sin2o; cos2z/ + cos2o; smzu ) + 

1 + e cos u 



2 e — sin z/ + sin^ i cos o; sin o; cos u + sin^ i cos^ u sin z/j ) 



sm z sin2o; (cos2z/ + e cosz/j ) 

1 + e cos u 



[l-e'^f^^ 2 A 2 •2- • o /"^^ cos 2z/ + e cos z^ 

n a V 1 — sm 2 sm 2o; / o- dz/ = (J . 



2 TT "'0 (1 + e cosz/) 

We see that the averages of all three Cartesian components of (^{df/d^l) x g — / x (9g/9i7) 
vanish. In the developments ( 11461 - IT^Sl) the following integrals were used: 



2'^ COSZ/ 1^ ~ ^ 3 TT e ri/io^ 

(1 + e cosz/)' " ~ V 1 + e (1 - e)' (1 + e)' ^ ^ 



and 



cos2z/ /I — e 3 vr e"' ,^ 

(150) 



(1 + e cosz/)' V 1 + e (1 - e)'' (1 + e) 
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A. 7. Calculation of /x(^xg-/x^ 

oi oi 



Just like the preceding subsections, this one is based on formulae ( IMll . ( l95ll . ( 1971) . and 
EED: 



5/ 



X 



df2 dh 



(151) 



1 + e cos V 



{ [ — sin(a; + v) sin 2 cos VL ] [ sin « cos {uj + v)] — 



cos z smiuj + y)] [ — sinfi sin {u + v) + cosQ cos(c(; + z/) cos i] } 



+ 



n a? a/1 - 
1 + e cosz/ 



e { [ — sin(u; + v) sin i cosfi ] [ sin i cosu ] — 



[ cos ^ sm{uj + u) ] [— sinfi sin a; + cosfi cosw cosz]} 



1 + e cos u 



I sin f2 sin^(co' + z/) cos z — cos f2 sin (a; + u) cos (a; + z/) 



+ 



e sin (a; + z/) [ sin fl sin cos i — cos cos u ] } 



dl 

di 



X 



dh dfi 
-rr- 9i ~ -TT- 93 



(152) 



1 + e cos z/ 



{ [ sin(co' + z/) cosz ] [— cosVL sin (cu + z^) — sinf2 cosz cos (a; + z/) ] — 
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[ sin Q sin i sin(a; + i^) ] [ sin i cos {uj + i/)] } 



+ 



n c? \J\ — 
1 + e cos nu 



e { [ sin(a; + i/) cosi] [— cosQ sin a; — sinQ cos a; cos z ] 



[ sin Q sin (a; + i/) sin z ] [ sin i cos a; ] } = 



1 + e cos 



— I — sin f2 sin(a; + v) cos(a; + i/) — cos f2 sin^ (a; + i/) cos i 



+ 



e [ — sin f2 sin (a; + v) cos io — cos ^2 cos i sin a; sin (a; + i/) ] } , 



df A dh dh 



oi 



(153) 



1 + e cos 



{ [ sin(a; + v) sin i sin Q ] [ — sin Q sin (a; + i/) + cos cos i cos(a; + v) ] 



[ — cos O sin(a; + i/) sin z ] [ — cos O sin (a; + i/) — sin O cos(a; + v) cos i ] } 



+ 



n a? \/l — 
1 + e COSI/ 



e { [ sin(a; + i/) sin i sin Q ] [ — sin Q sin ui + cos cos z cos ui ] 



[ — cos Q sin(a; + i/) sin i] [ — cos Q sin a; — sin Q cos a; cos i ] } = 



— r . o/ x • ■ • / \ • • • 1 

{ — sm (a; + I/) sm ^ — e smicu + i/) smo; sm z ^ , 

1 + e cos i/ ^ ^ ^ ' J 



Ol J ^ Ol Ol 



(154) 
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n 



1 + e COS!/ 



{ [ sin n cos(a; + z/) + cos fl sin(a; + u) cos i] cos(a; + u) cos i 



sm{u! + v) sin z [ — cos cos(a; + v) sin z ] } 



H e I [sinfi cosfo; + i/) + cosfi sinfo; + i/) cos il cos a; cos z 

1 + e cosi/ ^ ^ ^ ' ^ ' ^ 



sin(a; + z/) sin z [— cosfi cos a; sin z ] } = 



n 



1 + e cosv 



I cosVt sin (a; + u) cos (a; + i^) + sinil cos^(a; + u) cos z 



+ 



e [ cosr2 sin(6(; + v) cos a; + sinil cosiuj + i/) cos a; cos z ] } , 



J ^ — h -wr — Jl -wr 

Ol J ^ Ol Ol 



(155) 



n a 



1 + e cosu 



{ sin(a; + u) sin z sin ft sin z cos(a' + u) 



[ cos Q cos(a; + u) — sin Q sin(a; + u) cos z] cos(a; + u) cos z } 



+ 



n a 



1 + e cosz/ 



e { sin(a; + z/) sinz sinQ sin z coso; — 



[cosQ cos{uj + u) — sinQ sm{u! + u) cos z] coso^ cos z } 



n a" 



1 + e cos 



— I sinQ sm{u! + ly) cos{uj + v) — cosQ, cos^(a' + i/) cos z 
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+ 



e [sinf2 sm{u + u) cosu — cos 17 cos{uj + u) cosu cos z ] } , 



J X T7- — Jl -77- — J2 -Tr- 

ot J ^ 01 Ol 



(156) 



1 + e cos V 



{ [ cosfi cos(a; + z/) — sinf2 sin(c<; + v) cos i \ [ — cosfi sin i cos(a; + v) 



[ sinf2 cos((X' + z/) + cosf2 sin((X' + v) cos z] sinf2 cos((X' + v) cos i } 



+ 



nc? a/1 — 
1 + e cos z/ 



e { [cosfi cos(ci; + z/) — sinf2 sin(u; + z/) cos z ] [— cosfi sin i coso; 



sinf2 cos(a; + z^) + cosf2 sin(a; + z^) cos z] sinfi coso; cos « } 



n 



1 + e cosz/ 



I — sin z cos^(a; + z^) — e sin « cos(a; + z/) coso; | 



Now, by putting together, accordingly, fllSip with fll54p . fll52p with fllSSp . and fll53p with 
fll56p . we arrive, after some intermediate algebra, to the following expressions for the three 
Cartesian components: 



di 



(157) 



na'^y/l — 
1 + e cos u 



I sin cos i sin^(u; + z/) — cos^(w + z/) ^ — 2 cos sin (u + z/) cos {u + u) 



e [ sin Q cos i ( — cos {u + z/) cos + sin {u + z/) sin w ) — 2 cos sin {u + z/) cos w ] } 



1 + e cos u 



I ( — sin r2 cos « cos 2oj — cos f2 sin 2uj ) cos^ u — sin^ ^ ) 



44 



2 sin V cos v ( sin Q cos i sin 2ui 



— cos Q, cos 2a; ) ] + 



( — sin n cos i cos 2a; — cos O sin 2a; ) cos ^ + ( sin Vt cos i sin 2a; — 2 cos Vt cos^ a; j sin v j 



<9/ ^ r 

^ X g - / X — 



(158) 



^^^^^^^ — I [ cos n cos i i cos^fa; + i/) — sin^fa; + i/) ) — 2 sin Q sin iuj + v) cos (a; + i/) 

1 + e cos i/ L V ^ ^ ^ ^ ; ^ ' ^ ' 



+ 



e [ cos O cos i ( cos (a; + v) cos a; — sin (a; + z/) sin a; ) — 2 sin O sin (a; + v) cos lo ] } 



1 + e cos 1/ 



I 1^ ( cos Vt cos z cos 2a; — sin VL sin 2a; ) cos^ v — siv? ) + 



2 sin V cos ( — cos VL cos z sin 2a; — sin ^2 cos 2a; ) 1 + 



cos Vl cos i cos 2a; — sin Vt sin 2a; ) cos v + — 2 smVl cos^ a; + cos ^2 sin 2a; cos i ^ sin i/ 



9/ ^ -^g 

^ X g - / X — 



(159) 



V J- — c- f . . r 2 / \ • 2 / \ 

{ sm z cos [uj + u) — sm (a; + v) 

1 + e cos i/ L ^ ' ^ ^ 



+ 



e [ sin z cos a; cos(a; + z/) — sin z sin a; sin(a; + z/) ] } = 
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na^\/l — 
1 + e cos u 



I sin i cos 2u 



cos^ z/ — sin^ u 



2 sin 20; sini/ cosu 



+ 



e [ sin i cos 2u; cos z/ — sin i sin 2^ sin u]} . 

Averaging of the afore calculated three Cartesian components of (^{df/di) x g — / x (dg/di) 

almost exactly coincides with averaging of (^{df/d^l) x g — / x (dg/d^l)^ presented in 
the preceding subsection, and the result is the same, nil. Indeed, as can be easily seen from 
( 11571 - [TM]) . after we through out the (vanishing after averaging) inputs proportional to the 
odd functions sinu and sin2i/ , each of the three expressions (11571 - [T59|) will be propor- 
tional to (cos2i/ + e cosi/)/(l + e cosu) . Averaging of this expression via formula (I109p 
will vanish, as can be easily seen from (11491 - [T50|l . 



A. 8. Calculation of /x ( x g - / x 



dMo dMo 



This calculation will be the easiest: as 

df df dt ^ dt 



dMo dt dMo ° dMo 



and 



then, obviously. 



dMo dt dMo V l/r / ^^^o 



(160) 



dg dt ( Gm -\ dt 



A. 9. Calculation of ft- \ - f ^ ^ 

da 



With aid of f fTTOj) and ([68]) we establish: 

- fx — = - fx i (—] (—] 



- w A/Gma (1 - e2) - — = - w ^ ^ — — 

nl + ecosp\da) l + ecosv \da) 



(163) 
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where we used the formula 

f dt \ 1 — e cosE 1 1 — 



(164) 



\ dE / , , n n 1 + e cos u 

\ ' n, e, Mo 

that follows from the Kepler equation and from the first equation of (11051) . 
It remains to compute dE/da . From the Kepler equation 

E - esinE = Mo + n {t - to) (165) 

it follows that 

3 

dE (1 — e cos-E") = n a~'^ (t — to) da (166) 

wherefrom 

fdE\ _ 3 _^ n {t - to) _ 3 _^ E - e sinE - Mo 



da I , J., 2 1 — e cosE 2 1 — e cosE 

/ t, e, Mo 



(167) 



3 _i 1 + cos z/ 
- - n a- -Y^—r - to 



Insertion of the above into (11631) will give: 



/*x^ = ^wan(t-to) Vl - 



2 



aa 2 ^ ^ ^" ^ ' * " ' 

expression linear in t — to ■ 



A. 10. Calculation of fl - \ - f x ^ 



Repeating the line of reasoning presented in the preceding calculation, we derive from 
f lTT8l) . (I63D, and ffTMD the following: 

^ a/ ^ fdE\ ^ 2 (l-e^)'/' 

V^^/a,e,Af„ V^e;,^^^^,^ 1 + ecosz/ V^ey,^^^^^ 

To compute dE/de , one can start with the Kepler equation, 

E - e sin E = Mo + n {t - to) (169) 

which yields, for t, a, Mo being fixed: 

dE - e cos EdE - sinE de = , (170) 
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insertion of fllOSp wherein will entail 

/ dE \ sin E sin u 

\9^)t,a,M. ~ 1 - e COsE ~ VI - « 

All in all, 



(171) 



- / X ^ = - w sinz/ . . (172) 

de 1 + e cos v 



df A dh , df, 



(173) 



A. 11 Calculation of ii - \ - f ^ ^ 



Evidently, this vector is perpendicular to the instantaneous orbital plane and is uj- 
independent. A direct computation, based upon (p9]- 1101]) . will lead us to: 

Of A _ dh dh 

-TT- ^ J — TT- h — -TT- J 2 — 

OU j ^ OUJ OU 

[— sinf2 sin(a; + z/) + cosJl cos(a; + z/) cos 2 ] sin(ci; + z/) sin 2 — 
cos(a; + z/) sin z [sinfi cos(a; + z/) + cosi7 sin(a; + z/) cos z ] = 

a 2 [ — sm i 2 sm z J , 

(1 + e cosz/) 



5/ _ 9/3 9/i 

OU ] OU OU 

I 2 

(174) 



cos(a; + z/) sin i [cosfi cos(a; + z/) — sinf2 sin(ci; + z/) cos z ] — 
[— cosf2 sin(u; + z/) — sinf2 cos((X' + z/) cos i ] sin((X' + z/) sin « = 

2 (1 - e')' r n • ■ 1 
a 2 [ cos i Z sm z J , 

(1 + e cosz/) 



(175) 
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[— cosQ sm{uj + u) — sinQ cos{uj + u) cosi] [sinQ cos{u! + u) + cosQ sm{uj + u) cos i ] — 
[ — sin Q sin(a; + u) + cos Q cos(a; + i/) cos i ] [ cos Q cos{uj + u) — sin Q sin(a; + u) cos i ] 

— — r cos z = — a 2 cos z . 

(1 + e COSI/) 



Briefly, 



/x^] =-wa^^l^^ . (176) 
du (1 + e cosi/)' ^ ^ 



A.12. Calculation of fi-l -f x ^ 

oil 



on on ■'^ dn ■'^ 



[ cos Q cos(a' -\- u) — sin Q sin(a' + v) cos i ] sin(a; + v) sin z = 



(177) 



,2 



:i - e 



2\2 



1^1 + e COS!/) 



2 



[cosQ COS (a; + I/) — sinfi sin(a; + i/) cosi] sin(a; + i/) sin i , 



1 2 

(178) 

[ — sin O cos(a; + i/) — cos sin(a; + v) cos i ] sin(a; + v) sin z = 

(1 - e^f 

? 2 [sinQ cos(a; + I/) + cosQ sin(a; + i/) cosi] sin(a; + i/) sin i , 

(1 + e cosz/) 



(179) 
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[ — sin Q cos(a' + u) — cos Q sm{u! + u) cos i] [ sin Q cos(a' + u) + cos Q sin(a; + u) cos i ] — 
[cosQ cos{u! + u) — sinfi sin(a; + i/) cosi] [cosQ cos{uj + u) — sinQ sm{uj + u) cos i ] = 



— cos^(a' + u) — sin^(a; + u) cos^ i 



— a 



[1 + e cosu) 



cos^{uj + u) + sin^(a' + u) cos^ i 



A. 13. Calculation of fl i - f x 



df 

di 



dl 

di 



X f 



df2 r 

Ol 



Ol 



{ [ — cos Vt sin(a; + v) sin i ] sin(a' + v) sin i — 



sin(a; + u) cos i [ sin VL cos(a; + u) + cos sin(a; + v) cos i ] } 



(180) 



'1 - e 



2\2 



'1 + e cosz/) 



[— cosQ sm.{ijj-\-v) — sinQ cos(a; + i/) cos i ] sm.{uj -\- u) , 



df A _ df, df, 

Ol I Ol Ol 



{ sin(c<; + u) cos i [ cosVL cos{uj + u) — sinf2 sin(a; + u) cos i 



sinQ sin (a; + v) sin i sin(a; + v) sin i } 



(181) 



2^2 



(1 - e , 
(1 + e cosi/)' 



[— smVL sin (a; + I/) + cosO cos(a; + z/) cosz] sin(a; + i/) , 
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(182) 



[sinfi sin(a; + z/) sin z ] [sinfi cos(a; + z/) + cosfi sin(ci; + z/) cos z ] — 
[ — cosi7 sin((X' + z/) sin i ] [ cosi7 cos(u; + z/) — sinfi sin((X' + u) cos ^ ] 



(1 - e^f 

sinfw + v) cos(aj + z/) sin % = c? — 2 sin(cj + z/) cos(a; + v) sin z 



1 + e cos v) 



A. 14. Calculation of /7 I - / x "^-^ 



The Kepler equation (11681) entails that 

dj 1 df a 



3/2 



5Mo n dt y/G^ 
and, therefore, with aid of fl63l) we obtain: 



9/ 



7,3/2 



;i83) 



/ ^ 7^ = - ^ V^"^« (1 - = - w . (184) 



A. 15. The terms ( /x x / ) • ( /x x / ) 



In the course of numerical computation, it is convenient to keep such terms as parts of 
the Hamiltonian (l52l) . In case these terms are to be dealt with analytically, one will need 
the following general formulae valid for any Cj , j = a,e,fl,u,i, Mo : 



4(mx/>(,x/-).(M,/),(,,,-).L-,^ 



(/X X /) ■ /X 



X 



(185) 
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the neglect of the term with du/dCj being legitimate, because this input is of a even 
higher order than the other one 1^1 A direct algebraic calculation then yields, under the said 
approximation: 



expression for fj given above by ([99] - IIOII) 

The physical meaning of dft/dCj is the influence of the satellite upon the precession rate of the 
primary. As demonstrated by Laskar (2004), this effect is of a higher order of smallness. 



52 



References 



[I] Arnold, V. I. 1989. Mathematical Methods of Classical Mechanics. Springer- Verlag, NY 
1989. 

[2] Ashby, N., and T. Allison. 1993. "Canonical Planetary Perturbation Equations for 
Velocity-Dependent Forces, and the Lense-Tliirring Precession." Celestial Mechanics and 
Dynamical Astronomy, Vol. 57, pp. 537 - 585. 

[3] Brouwer, D. 1959. "Solution of the problem of artificial satellite theory without drag." 
The Astronomical Journal, Vol. 64, pp. 378 - 397. 

[4] Brouwer, D., and G. M. Clemence. 1961. Methods of Celestial Mechanics. Academic 
Press, NY & London. (Chapter XI) 

[5] Brumberg, V. A., L. S. Evdokimova, & N. G. Kochina. 1971. "Analytical Methods for 
the Orbits of Artificial Satellites of the Moon." Celestial Mechanics, Vol. 3, pp. 197 - 221. 

[6] Brumberg, V.A. 1992. Essential Relativistic Celestial Mechanics. Adam Hilger, Bristol. 

[7] Burns, J., & V.Safronov 1973, "Asteroid Nutation Angles." The Monthly Notices of the 
Royal Astronomical Society of London, Vol. 165, pp. 403 - 411 

[8] Chambers, J. E. 1999. "A Hybrid Symplectic integrator that Permits Close Encoun- 
ters between Massive Bodies." The Monthly Notices of the Royal Astronomical Society of 
London, Vol. 304, pp. 793 - 799. 

[9] Chernoivan, V. A., & Mamaev, I. S. 1999. "The Restricted Two-Body Problem and the 
Kepler Problem in the Constant-Curvature Spaces." Regular and Chaotic Dynamics, Vol. 
4, pp. 112 - 124. 

[10] Defraigne, P.; Rivoldini, A.; Van Hoolst, T.; Dehant, V. 2003. "Mars nutation resonance 
due to Free Inner Core Nutation." Journal of Geophysical Research, Vol. 108, p. 5128 

[II] Dehant, V.; van Hoolst, T.; Defraigne, P. 2000. "Comparison between the nutations of 
the planet Mars and the nutations of the Earth." Surveys in Ceophysics, Vol. 21, pp. 89 
- 110. 

[12] Duncan, M. J.; Levison, H. F. & M. H. Lee. 1995. "A Muhiple Time Step Symplectic 
Algorithm for Integrating Close Encounters." The Astronomical Journal, Vol. 116, pp. 
2067 - 2077 

[13] Efroimsky, M. 2001. "Relaxation of Wobbling Asteroids and Comets. Theoretical Prob- 
lems. Perspectives of Experimental Observation." Planetary & Space Science, Vol. 49, 
pp. 937 - 955 

[14] Efroimsky, M. 2002a. "Equations for the orbital elements. Hidden symmetry." Preprint 
No 1844 of the Institute of Mathematics and its Applications, University of Minnesota 
|http: / / www.ima.umn.edu/preprints/ feb02 / feb02.html| 

[15] Efroimsky, M. 2002b. "The Implicit Gauge Symmetry Emerging in the N-body Problem 
of Celestial Mechanics." astro-ph/02 12245. 



53 



[16] Efroimsky, M. 2002c. "Euler, Jacobi, and Missions to Comets and Asteroids." Advances 
in Space Research, Vol. 29, pp. 725 - 734 

[17] Efroimsky, M., & P. Goldreich. 2003. "Gauge SjTiimetry of the N-body Problem in the 
Hamilton- Jacobi Approach." Journal of Mathematical Physics, Vol. 44, pp. 5958 - 5977 
|astro-ph/0305344| 

[18] Efroimsky, M., & Goldreich, P. 2004. "Gauge Freedom in the N-body Problem of Ce- 
lestial Mechanics." Astronomy & Astrophysics, Vol. 415, pp. 1187 - 1199 
|astro-ph/0307130, 

[19] Eubanks, T. M. 1993. "Variation in the orientation of the Earth. In: Contributions 
of Space Geodesy to Geodynamics: Earth Dynamics, pp. 1-54, Amer. Geophys. Union, 
Washington, 1993. 

[20] Goldreich, P. 1965. "Inclination of satellite orbits about an oblate precessing planet." 
The Astronomical Journal, Vol. 70, p. 5 - 9 

[21] Gradshtein, I. S., and I. M. Ryzhik. 1994. Tables of Integrals, Series, and Products. 
Academic Press, NY. 

[22] Hill, G. W. 1913. "Motion of a System of Material Points under the Action of Gravita- 
tion." The Astronomical Journal, Vol. 27, p. 171 - 182 

[23] Kinoshita, T. 1993. "Motion of the Orbital Plane of a Satellite due to a Secular Change 
of the Obliquity of its Mother Planet." Gelestial Mechanics and Dynamical Astronomy, 
Vol. 57, pp. 359 - 368 

[24] Kozai, Y. 1960. "Effect of precession and nutation on the orbital elements of a close 
earth satellite." The Astronomical Journal, Vol. 65, pp. 621 - 623. 

[25] Lagrange, J.-L. 1778. Sur le Probleme de la determination des orbites des cometes 
d'apres trois observations, 1-er et 2-ieme memoires. Nouveaux Memoires de I'Academie 
de Berlin, 1778. Later edition: In CEuvres de Lagrange. Vol. IV, Gauthier-Villars, Paris 
1869. 

[26] Lagrange, J.-L. 1783. Sur le Probleme de la determination des orbites des cometes 
d'apres trois observations, 3-ieme memoire. Ibidem, 1783. Later edition: In Q^uvres de 
Lagrange. Vol. IV, Gauthier-Villars, Paris 1869. 

[27] Lagrange, J.-L. 1808. "Sur la theorie des variations des elements des planetes et en 
particulier des variations des grands axes de leurs orbites" . Lu, le 22 aout 1808 a I'lnstitut 
de France. Later edition: In QSuvres de Lagrange. Vol. VI, p. 713 - 768, Gauthier-Villars, 
Paris 1877. 

[28] Lagrange, J.-L. 1809. "Sur la theorie generale de la variation des constantes arbitraires 
dans tons les problemes de la mecanique". Lu, le 13 mars 1809 a I'lnstitut de France. 
Later edition: In CEuvres de Lagrange. Vol. VI, p. 771-805, Gauthier-Villars, Paris 1877. 

[29] Lagrange, J.-L. 1810. "Second memoire sur la theorie generale de la variation des con- 
stantes arbitraires dans tons les problemes de la mecanique". Lu, le 19 fevrier 1810 a 
I'lnstitut de France. Later edition: In CEuvres de Lagrange. Vol. VI, p. 809 - 816, Gauthier- 
Villars, Paris 1877. 



54 



[30] Laskar, J. 1990. "The Chaotic Motion of the Solar System. A Numerical Estimate of 
the Size of the Chaotic Zones." Icarus, Vol. 88, pp. 266 - 291. 

[31] Laskar, J. 2004. "A Comment on "Accurate Spin Axes and Solar System Dynamics: 
Climatic Variations for the Earth and Mars"." Astronomy & Astrophysics, Vol. 416, pp. 
799 - 800. 

[32] Laskar, J., & J. Robutel. 1993. "The Chaotic Obhquity of the Planets." Nature, Vol. 
361, pp. 608 - 612. 

[33] Marsden, J., & T. Ratiu. 2003. Introduction to Mechanics and Symmetry. Springer, NY 
2003 

[34] Molina, A., F. Moreno, & F. Martinez-Lopez. 2003. "Energy dissipation by internal 
stresses in a free-rotating symmetric ellipsoid: Application to Comet P/Halley." Astron- 
omy & Astrophysics, Vol. 398, pp. 809 - 817 

[35] Morbidelli, A. 2002. Modem Celestial Mechanics: Dynamics in the Solar System. Taylor 
& Francis, London 2002. 

[36] Murison, M. 1988. "Satellite Capture and the Restricted Three-Body Problem." Ph.D. 
Thesis, University of Wisconsin, Madison 1988. 

[37] Newman, W., & M. Efroimsky. 2003. "The Method of Variation of Constants and Mul- 
tiple Time Scales in Orbital Mechanics." Chaos, Vol. 13, pp. 476 - 485. 

[38] Pang, K. D., J. B. Pollack, J. Vcvcrka, A. L. Lane, & J. M. Ajello. 1978. "The Compo- 
sition of Phobos: Evidence for Carbonateous Chondrite Surface from Spectral Analysis." 
Science, Vol. 199, p. 64 

[39] Pollack, J. B.; Burns, J. A.; & Tauber, M. E. 1979. "Gas Drag in Primordial Circum- 
planetary Envelopes. A Mechanism for Satellite Capture." Icarus, Vol. 37, p. 587 

[40] Prendergast, K. H. 1958. "The Effects of Imperfect Elasticity in Problems of Celestial 
Mechanics." Astronomical Journal, Vol. 63, pp. 412- 414 

[41] Proskurin, V. F., & Batrakov, Y. V. 1960. "Perturbations of the Motion of Artificial 
Satellites, caused by the Earth Oblateness." The Bulletin of the Institute of Theoretical 
Astronomy, Vol. 7, pp. 537 - 548. 

[42] Richardson, D. L., and Kelly, T. J. 1988. "Two-body motion in the post- Newtonian 
approximation." Celestial Mechanics, Vol. 43, pp. 193-210. 

[43] Sharaf, S. G. & N. A. Budnikova. 1969a. "On Secular Perturbations in the Elements of 
the Earth's Orbit and their Influence on the Climates in the Geological Past." Bulletin of 
the Leningrad Institute of Theoretical Astronomy. Vol. 11, pp. 231 - 261. /in Russian/ 

[44] Sharaf, S. G. & Budnikova. 1969b. "Secular Perturbations in the Elements of the Earth's 
Orbit and the Astronomical Theory of Climate Variations" Proceedings of the Leningrad 
Institute of Theoretical Astronomy. Vol. 14, pp. 48 - 84. /in Russian/ 



55 



[45] Sharma, I., J. A. Burns, & C.-Y. Hui. 2004. "Nutational Damping Times in Solids of 
Revolution." To be published in the Monthly Notices of the Royal Astronomical Society 
of London. 

[46] Slabinski, V. 2003. "Satellite Orbit Plane Perturbations Using an Efroimsky Gauge 
Velocity." Talk at the 34th Meeting of the AAS Division on Dynamical Astronomy, Cornell 
University, May 2003. 

[47] Szebehely, V. 1967. Theory of Orbits. Academic Press, NY. 

[48] Tolson, R. H., and 15 collaborators. 1978. "Viking First Encounter of Phobos. Prelimi- 
nary Results." Science, Vol. 199, p. 61 

[49] Touma, J., & J. Wisdom. 1994. "Lie-Poisson Integrators for Rigid Body Dynamics in 
the Solar System." The Astronomical Journal, Vol. 107, pp. 1189 - 1202. 

[50] Van Hoolst, T.; Dehant, V.; Defraigne, P. 2000. "Chandler wobble and Free Core Nu- 
tation for Mars." Planetary and Space Science, Vol. 48, pp. 1145 - 1151. 

[51] Vernekar, A. D. 1972. "Long-period Global Variations of the Incoming Solar Radiation." 
Meteor. Monogr., Vol. 12, pp. 1 - 21. 

[52] Veverka, J. 1977. "Phobos and Deimos." Scientific American, Vol. 236, p. 30 

[53] Ward, W. 1973. "Large-scale Variations in the Obhquity of Mars." Science. Vol. 181, 
pp. 260 - 262. 

[54] Ward, W. 1974. "Climatic Variations of Mars. Astronomical Theory of Insolation." 
Journal of Geophys. Research, Vol. 79, pp. 3375 - 3386. 

[55] Ward, W. 1979. "Present Obliquity Oscillations of Mars - Fourth-order Accuracy in 
Orbital e and i." Journal of Geophys. Research, Vol. 84, pp. 237 - 241. 

[56] Ward, W. 1982. "Comments on the long-term Stability of the Earth's Obliquity." Icarus, 
Vol. 50, pp. 444 - 448. 

[57] Waz, P. 1999. "Disturbing Function in the Analytical Theory of the Motion of Phobos." 
The Astronomical Journal, Vol. 348, pp. 300 - 310 



56 



